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INTRODUCTION 


This research project represented a cooperative 
effort between Lewis Research Center (LeRC) and Cleveland 
State University (CSU). This project has been continued 
under contract with Analex Corporation. 

SUMMARY OF RESEARCH ACCOMPLISHMENTS 

Initial investigation by the principal investigator 
occurred during the summer of 1988 under a NASA/ASEE Summer 
Faculty Fellowship, and discussed in the Interim Status 
Report *1. Continued investigation which occurred under 
NAG 3 - 1 0 0 8 , for the time period February 8, 1989-June 15, 
1989, was also reported in the first interim status report. 

Additional work was continued during Summer, 1989, 
under a NASA/ASEE Summer Faculty Fellowship. This was also 
Included in the first interim status report. 

NASA LeRC approved the continued funding of this 
grant for the period September 15, 1989 through September 
14, 1990. The following accomplishments occurred during the 
period September 15, 1989 through March 15, 1990: 
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1. A rigorous solution for the dynamic analysis of a 
free/free beam with an axial tension pre-load was 
developed. Solution of the characteristic equation 
suggested that the three required rigid body modes were 
present . 

2. A paper en t i t 1 ed , " Dynam i c Analysis of Space-related 
Linear and Non-linear Structures", was co-authored by 
Professor Bosela, Dr. Francis Shaker (NASA LeRC), and 
Professor Demeter Fertis (University of Akron, Akron, 
Ohio). The paper was presented by Professor Bosela at 
the Southeast Conference of Theoretical and Applied 
Mechanics XV. Atlanta, Georgia, during March, 1990. 

3. A three-node beam element was developed, using Martin's 
methodology. It was duplicated using a variational 
formulation. Its performance in various sample problems 
was tested. 

4. The bow-string problem was identified as an idealized 
model of a solar array which had potential to yield 
rigid body rotation capability. 

5. Numerous papers regarding derivation of higher-order 
stiffness matrices were reviewed, including work by 
Argyris, Saunders, Paz , Martin, Marcal, and others. 

Their matrices were tested for rigid body capabilities. 
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From March 16. 1990, through the completion of this 


grant, the following accomplishments were made: 


1. Exact solutions of various pre-loaded beam problems were 
examined, and the Galerkin criterion was used to develop 
stiffness and mass matrices. 

2. The modified matrices developed using the Galerkin 
criterion were incorporated into a finite element 
dynamic analysis algorithm, and the resulting finite 
element solution compared with the rigorous solution. 

3. A directed force correction matrix for the pre-loaded 2 
dimensional beam element was developed at the global 
level. This matrix produced a tangential stiffness 
matrix which does possess all of the required rigid body 
modes . 

4. This global force correction was incorporated into a 

finite element dynamics algorithm, and was shown to 
correct the missing zero eigenvalue customary in 
traditional finite element solutions, without affecting 
the eigenvalues corresponding to the flexible modes. It 
also performed very well in the 

di agona 1 i zat i on/par t i t i on i ng methodology used in matrix 
dynamic analysis. 
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5. The detailed results of this study were published as l 
doctoral dissertation at the University of Akron. Akron, 
Ohio, and are attached. The period of performance for 
this grant expired on April 6, 1991. The final 

dissertation was submitted for publication during the 
summer of 1991. and is available at the University of 
Akron library. 
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ABSTRACT 


Space structures, such as the space station solar 
arrays, must be extremely light-weight, flexible structures. 
Accurate prediction of the natural frequencies and mode 
shapes is essential for determining the structural adequacy 
of components, and designing a controls system. The tension 
preload in the "blanket" of photovoltaic solar collectors, 
and the free/free boundary conditions of a structure in 
space, causes serious reservations on the use of standard 
finite element techniques of solution. In particular, a 
phenomena known as "grounding", or false stiffening, of the 
stiffness matrix occurs during rigid body rotation. 

This dissertation examines the grounding phenomena 
in detail. Numerous stiffness matrices developed by others 
are examined for rigid body rotation capability, and found 
lacking. Various techniques are utilized for developing new 
stiffness matrices from the rigorous solutions of the 
differential equations, including the solution of the 
directed force problem. A new directed force stiffness 
matrix developed by the author provides all the rigid body 
capabilities for the beam in space. 
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CHAPTER 1 


INTRODUCTION 

In order to be cost-effective, space structures must 
be extremely light-weight, and subsequently, very flexible 
structures. The power system for Space Station Freedom is 
such a structure. Each array consists of a deployable truss 
mast and a split "blanket" of photovoltaic solar collectors. 
The solar arrays are deployed in orbit, and the blanket is 
stretched into position as the mast is extended during 
deployment. Geometric stiffness due to the tension preload 
in the blanket make this an interesting non-linear problem. 

The space station will be subjected to various 
dynamic loads, during shuttle docking, solar tracking, 
attitude adjustment, etc.. Accurate prediction of the 
natural frequencies and mode shapes of the space station 
components, including the solar arrays, is critical for 
determining the structural adequacy of the components, and 
for designing a dynamic controls system. 

This dissertation has the following objectives: 

1. Examine in detail the "grounding" phenomenon associated 
with rigid body rotation of a pre-loaded beam in space. 
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2. Examine beam geometric stiffness matrices developed by 
others with respect to rigid body motion capabilities. 

3. Develop higher order stiffness matrices from the 
rigorous solution utilizing Galerkin's criterion, 
incorporate these stiffness matrices into a finite 
element algorithm, and compare the finite element 
solutions with the rigorous solutions. 

4. Examine the directed force (bow-string) problem for its 
potential as a basis for developing stiffness matrices 
which possess rigid body rotational capabilities. 

5. Check the performance of any new matrix which possesses 
a complete set of rigid body motion capabilities in the 
diagonal ization/partitioning methodology used in dynamic 
response . 


CHAPTER 2 


LIMITATIONS OF CURRENT METHODOLOGY 

Most structural systems are rigidly attached to 
supports at either or both ends. In order for any movement 
to occur, the structure must deform, and internal strain 
energy is developed. Space structures, on the other hand, 
are not rigidly attached to the ground. Instead, they are 
free to move as rigid bodies as well as to deform. 

Complex structures are generally analyzed using 
finite element computer programs which solve the dynamic 
equations of motion using matrix analysis techniques. The 
equations of motion are set up in the form of the 
generalized eigenvalue problem 

[[K]-Qi 2 [M]] {Ui> = {Ri) 

where [K] is the global stiffness matrix 
[M] is the global mass matrix 
Qi are the natural frequencies of vibration 
{ Ui } are the displacement or mode shape vectors 
{R^} are the forces 
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Using that basis, rigid body inodes are the eigenvectors 
associated with zero frequencies of vibrations 
(eigenvalues) . 

Current methodology utilizes MSC/NASTRAN solution 64 
to generate the tangential stiffness matrix for the deployed 
array, storing this matrix in a database, then using this 
matrix in solution 63 dynamic analysis, to obtain the 
frequencies of vibration. As a routine check of the model, 
the global stiffness matrix is multiplied with a matrix of 
the rigid body modes to determine whether any pseudo-forces 
occur. (Whether strain energy has developed.) Since no 
internal stresses should occur during rigid body motion, the 
generation of pseudo-forces indicates that an internal 
"grounding", or false stiffening, of the system occurs, due 
to errors or deficiencies in the finite element model. 

It was found that the global stiffness matrix does 
not possess rigid body rotation capabilities. In order to 
predict the dynamic response of the structure, a Craig- 
Bampton substructuring scheme is used. However, certain 
erroneous non-zero terms appear in the null set of the 
partitioned matrices due to the grounding effect. They must 
be zeroed out, and the missing rigid body modes appended to 
the matrix, in order to more accurately predict the dynamic 
response [1]. 


5 


The author idealized the problem as a free/free beam 
in tension, and found [2] that the pseudo-forces are 
developed at the element level due to limitations inherent 
in the geometric stiffness matrices currently in acceptable 
use. In particular, the geometric stiffness matrices for 
the beam element lack the capability for rigid body 
rotations, especially when the rotations are large. 

The geometric (initial stress) stiffness matrices in 
current use developed from a Bernoul 1 i-Eul er formulation 
have been shown to provide acceptable results for most 
static displacement and buckling problems, provided a 
sufficient number of elements are used [3]. However, 
refinement of the mesh does not produce convergence to the 
missing zero frequency in the dynamics problem of the pre- 
loaded beam with free/free boundary conditions. In 
addition, higher frequencies may be significantly in error. 
Table 1 compares the finite element solution for a pre- 
tensioned beam with pinned/rol ler and free/free boundary 


conditions . 
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TABLE 1 Comparison of Finite Element Method Versus Exact 
Solution for a Beam in Tension 

A = 48 in^ 

E = 30 x 10^ psi 
I = 1000 in 4 

2 2 

m = 0.03525 lb-sec /in 
P = 10,000,000 lb 
L = 100 IN 


f req 

number of 
1 2 

elements 

4 

8 

% error 

rigorous sol 

. 19 
Fertis 

■ 

1142 

1056 

1053 

1053 

3.5 

1043 


3501 

3257 

3195 

3180 

axial 


H 

4758 

4180 

3807 

3794 

4.0 

3647 

■ 


10291 

8494 

8357 

9.0 

7669 





Free/Free 

Beam 


f req 

number of 
1 2 

elements 

4 

8 

% error 

rigorous sol 

16 17 

Bosela/Shaker 

■ 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


580 

580 

579 

579 

to 

0 

B 

2798 

2383 

2381 

2378 

15.0 

2017 

B 

7001 

6725 

5991 

5958 

4.7 

5725 


2 

3 


Pin/Rol 1 er 












CHAPTER 3 


ELASTIC STIFFNESS MATRIX 


The elastic stiffness matrix for a 2-node Bernoulli 

beam is 


[Ke] = El 


AL 2 /I 


-al 2 /i 


0 

12 

6L 

0 

•12 


6L 

4L 2 

0 

-6L 


-al 2 /i 


al 2 /i 


0 

-12 
-6L 2L 
0 


6L 

2 


12 -6L 


6L 


2L‘ 


-6L 4L‘ 


The [Ke] matrix must possess the capacity of a full 
set of rigid body modes. In other words, the element must 
be able to both translate and rotate without developing 
stresses (see Figure 1). 
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Rigid Body Rotation 


Figure 1 


c . 

{^RBR} = 0 [0 , “L/2 , 1 ,0 ,L/2 , 1] 


Rigid Body Modes 
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Note that in Fig 1(c) that the rotation is 
considered to be relatively small, such that the 
displacement in the axial direction due to the rotation is 
negl igibl e . 

Multiplying [Ke] x [Rigid Body Mode] = 

[ 0 , 0 , 0 , 0 , 0 , 0 ] 

holds for all three modes. Hence, [Ke] possesses all the 
required rigid body mode capabilities. 

Another way of determining whether [Ke] possesses 
all the rigid body mode capabilities is to solve the 
dynamic analysis of the beam with free/free boundary 
conditions. This was done [5] using the finite element 
dynamics algorithm in the computer program NLFINITE . FOR 
(Appendix C) . The results were three zero eigenvalues and 
corresponding rigid body mode shapes. 

Another beam stiffness matrix which incorporates 
shear effects is referred to as a Timoshenko beam. The 
elastic stiffness matrix for a Timoshenko beam is 


EI(l/(l+5)) 
[Ke] = = 


AL 2 ( 1 + 5 )/ I 0 
0 12 

0 6L 

-AL 2 ( 1 + 5 ) / I 0 
0 -12 

0 6L 


0 -AL 2 (l+5)/I 0 0 

6L 0 -12 6L 

( 4+5 )L 2 0 -6L (2-5)L 2 

0 AL 2 ( 1+5 ) / I 0 0 

-6L 0 12 -6L 

( 2-5)L 2 0 -6L ( 4+5 )L 2 
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Where $ = 12 EI/(L 2 K'AG), which corrects for shear 
deformation. As K'AG becomes very large, 1*0, and [Ke]<y = 
[Ke] . The Timoshenko elastic stiffness also possesses a 
full set of rigid body modes. 

A major difference in the Timoshenko approach is 
that the bending rotation is considered independently in the 
derivation, not simply the derivative of the displacement 
equation, as is done in the Bernoulli derivation. 


CHAPTER 4 


GEOMETRIC STIFFNESS MATRIX DEVELOPMENT 

The presence of an axial force introduces additional 
stiffness terms, resulting in the geometric stiffness, or 
initial stress stiffness matrix. Various formulations of 
the geometric stiffness matrix have been developed. 

When the Hermitian interpolating polynomials (used 
to derive the [Ke] matrix) are used in deriving the 
geometric stiffness coefficients [4], the resulting [Kg] is 
referred to as the consistent geometric stiffness matrix 
(Bernoulli beam geometric stiffness). 

0 0 0 0 0 

36 3L 0 -36 3L 

3L 4L 2 0 -3L -L 2 

0 0 0 0 0 

36 -3L 0 36 -3L 

3L -L 2 0 -3L 4L 2 

Application of the rigid body modes to [Kg] results in 
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12 


[Kg] x 


1 

0 

0 


0 

0 

0 

0 

1 

-L0/2 


0 

0 

-P8 

0 

0 

e 

= 

0 

0 

0 

1 

0 

0 


0 

0 

0 

0 

1 

L0/2 


0 

0 

pe 

0 

0 

0 


0 

0 

0 


The terms ± P0 are fictitious forces generated 
during the rigid body rotation. Similarly, dynamic 
analysis, using NLF I N I TE . FOR , yields only two zero 
eigenvalues for the free/free beam in tension, corresponding 
to axial and transverse rigid body translations only. 

Various formulations have been used for establishing 
the geometric stiffness matrices from the static 
displacement problem. Martin [6] used a strain energy 
formulation with interpolating polynomials. Clough [7] used 
minimization of the potential function with the Hermitian 
polynomials. Both approaches yield a consistent geometric 
stiffness matrix, which lacks rigid body rotation 
capability, as was previously demonstrated. 


The author followed Martin's methodology in 
developing a 3-node beam geometric stiffness matrix [8]. 
The following matrix was obtained. 
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3L 


n 

18272 







u 

n 

105L 

3469 

659L 




SYMMETRIC 

u 

105 

105 






-20 

3L 

0 

0 

64 

3L 





n 

-22096 

-4208 

n 

27292 




u 

105L 

105 

u 

105L 




0 

304 

23L 

A 

-72 

109L 



5 

2 

U 


5 



13 

3L 

0 

0 

-44 

3L 

0 

0 

31 

3L 


0 

3824 

739 


-5296 

56 

n 

1472 

105L 

105 

U 

105L 

5 


105L 

0 

-3469 

-659L 

0 

4208 

-23L 

fi 

-739 659L 

105 

105 

105 

2 


105 105 


[Kg] 

3 -NODE 

By inspection, [Kg] 3 _ no( j e has two rigid body 
translational capabilities. The exact rigid body rotation 
vector is 

[L(l-COS(2B))/2,-LSIN(2B)/2,2B,0,0,2fl, 

-L ( l-COS( 2fl))/2 / LSIN(2B)/2,2B] T , 
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where B is 1/2 the angle of rotation. 

If this vector is expanded in power series form, 
upon retaining the first two terms, and factoring out BL, 
one obtains 

{ u RBR> T = [fi-B 3 /3,-l+2fi 2 /3,2/L,0,0,2/L,-fi+B 3 /3,l-2B 2 /3,2/L] . 

Multiplying [Kg] by {Urjjr} yields 

[-2fl 2 ,91.73B 3 -16B, l(17.3fl 3 -3B) , 8fl 2 , 16B-106 . 7B 3 , 

L(33.01B 3 -6B),-6B 2 ,14.93B 3 ,L(3B-17.33B 3 )], 

which contains numerous non-zero terms. Hence, [Kg] 3 _ no( i e 
does not possess rigid body rotation capability. 

Saunders [9] solves for the exact solution of the 
differential equations using a Timoshenko approach, then 
expands his "exact" stiffness matrix in a power series 
solution, obtaining a series of matrices of increasing 
order . 

Saunders "exact" stiffness matrix is 


P 

[K] =— 
Z 


BR-SIN(BL) 

SYMMETRIC 

SIN(BL) 

1-COS BL -L-COS(BL) 

BR 

-BR-SIN(BL) COS ( BL ) -1 BR-SIN(BL) 

SIN(BL) SIN(BL) 

l-COS(BL) L COS(BL)-l L-COS(BL) 


BR 


BR 
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where 

B = V ( P/EIR) 
a = BL 

R = (l-P/K'AG) 

P = AXIAL LOAD 

K'AG = beam shear rigidity 

I = moment of inertia 

z = SIN(BL) • (2* TAN ( BL/ 2 ) -BLR) 

By observation, rigid body translation capability is 
present in the transverse direction. 


Upon multiplying [K]*{Urbr) = 


P0 


- 2COS ( a ) - LRB *SIN(a)+2 
0 

2COS(a)+LRB-SIN(a)-2 

0 


For small a, COS(a)s:l, SIN(a)s:ct. 
-2COS(a)-LRB-SIN(a)+2 = -LRfia 

= -LRB • BL 
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P 0 P0 

z sin(a) (2'Tan (a/2) - aR) 


P0 

a (a - aR ) 


P0 

a 2 (1 - R) 


Thus, (P/z) (-PL 2 /EI) = P0(-B 2 L 2 R) 

a 2 (1-R) 

- P0R 
(1-R) 

- P 0 ( K ' AG - P) 

K 1 AG 
P 

K'AG 

= (P - K'AG)0. 


Thus, Saunders' "exact stiffness matrix does not 


possess the required rigid body rotation mode. 
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Argyris [10] uses his "natural formulation" to 
develop [Ke] and [Kg], which are identical with traditional 
[Ke] and [Kg]. He obtains another matrix [K nc ], referred to 
as his load correction matrix, which compensates for non- 
conservative forces. 

If we consider the axial load to remain tangent to 
the slope of the beam at the end points, Argyris's total 
geometric stiffness matrix [Kg + K nc ] becomes 



0 0 

SIN2B 

0 

0 

0 


0 6/5L 

1/10+COS2B 

0 

-6/5L 

1/10 

[Kg] = P 

0 1/10 

2L/15 

0 

-1/10 

-L/30 

TOTAL 

0 0 

0 

0 

0 

-SIN2B 


0 -6/5L 

-1/10 

0 

6/5L 

-1/10-COS2B 


0 1/10 

-L/30 

0 

-1/10 

2L/15 

The matrix is : 

nonsymmetric . 


Multiplying [Kg ] total 

by the exact 

rigid body 

rotation vector, then 

applying small 

angle considerations. 

2 

yields [4B 

,0 

, 0 , - 4B 2 

T 

,0,0] , which 


contains non-zero terms. Hence, [K9 ]t 0TAL does not possess 
rigid body rotation capability. Note that the pseudo-forces 
now occur in the axial direction. 

Martin [11] summarizes work done by Marcal [12] 
which introduced higher order terms in his initial 
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displacement matrices. In addition to the conventional [Ke] 
and [Kg], his initial displacement matrices are 

0 0 0 ”^4 0 

b4 b2 0 "^*4 ~b2 0 

0 0 0 0 0 0 

[»!]=AE/L 

0 “^4 0 0 b 4 0 

-b4 -b2 0 b4 b2 0 

0 0 0 0 0 0 


and 


0 0 0 0 

0 0 -1 . 5b 4 2 0 

0 0 0 0 

0 0 0 0 

0 0 1 . 5b 4 2 0 

0 0 0 0 

where u = b^ + b 2 X 
v = b 3 + b 4 X. 
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The basic non-linear equation is 

[K + 1/2 + 1/3 fl 2 ]*{U} = (R>- 

By inspection, and [ S 2 3 possess the required 

rigid body translation capabilities. 

Let the rotation angle = 2B (Figure 2). 

ui = hi vi=b 3 

U2=bi + b2L V2 = b 3 + b 4 L 

b 4 = (v 2 -vi)/L 
b 4 = (LB + LB) /L 
b 4 = 2fl 

b 4 2 /2 = 4B 2 /2 = 2fl 2 . 

2 

Similarly, b 2 = -2B . 

Thus, the rigid body rotation check becomes 
[K + + N 2 ] - {Urbr} must equal 0, 

where 

B 0 0 -B 0 

-B 2 0 -B B 2 0 

0 0 0 0 0 

-B 0 0 B 0 

B 2 0 B -B 2 0 

0 0 0 0 0 


0 

B 

0 

[N X ] =AE/L 

0 

-B 

0 
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and 

0 0 0 0 0 
2B 2 0 0 -2B 2 0 

0 0 0 0 0 

0 0 0 0 0 

2B 2 0 0 2B 2 0 

0 0 0 0 0 

Performing the rigid body rotation check yields 

[ 0 , 4AB 3 E- 4B 3 LP-2BP ,0,0, -4AB 3 E+4B 3 LP+2BP , 0] T . 

Note that non-zero pseudo-force terms still appear. 

Development of the stiffness matrices from the 
equation of motion has been investigated by Paz, using both 
a Bernoulli [13] and Timoshenko [14] beam approach. He 
developed his "exact" stiffness matrix, then expanded it in 
a power series solution. 

His solution, based on the transverse vibration of 
a beam with an axial compression load, is of the form 

[s] = [K] - [G 0 ] P - [M 0 ]Q 2 - [A^PQ 2 - [GqJP 2 - [M^ ] Q 4 •••. 



where 
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[K] is the traditional elastic stiffness matrix with no 
axial terms. 

[Gq] is the standard geometric stiffness matrix. 

[Mq] is the first order mass matrix (consistent mass 
matrix) . 


[Mq ] = mL/ 420 


[A^] is the 


[A 1 ]=mL 3 /EI 


[G]J is the 


156 

symmetric 


22L 

4L 2 


54 

13L 156 


-13L 

-3L 2 -22L 4L 2 


second o 

rder mass-geometrical matrix 

1/3150 

SYMMETRIC 

L/1260 

L 2 /3150 


-1/3150 

L/1680 1/3150 

-L/1680 

L 2 / 3600 -L/1260 L 2 /3150 

second order geometrical matrix. 


[Gi] = 1/EI 


L/700 

L 2 /1400 

-L/700 

L 2 /1400 


. SYMMETRIC 

11L / 6300 

-L 2 /1400 L/700 

-13L 3 /12600 -L 2 /1400 llL 3 /6300 
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[H^] is the second order mass matrix. 



59 


SYMMETRIC 



161-7 





223L 

71L 2 



2 t 5 
m L 

2910-6 

1279 

4365-9 

1681L 

59 


1000 El 


3880-8 

23284-8 

161-7 



-1681L 

-1097L 2 

-223L 

71L 2 


2384-8 

69854-4 

2910-6 

4365-9 


The mass matrices don't possess rigid body modes, 
but they are not intended to, since they generate the 
inertial forces. [G^] possesses all the rigid body modes. 
Hence, no correction to [Go], which lacks rigid body 
rotation capability, is applied. Thus, "grounding" during 
rigid body rotation still occurs. 

Similarly, Paz's Timoshenko formulation (which 
includes rotary inertia and shear terms), generates the 
matrix 


[R 0 ] =mL/30 ) (R/L) 2 ( 1+E/K ' G) 
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SYMMETRIC 

3L 

CM 


-36 

-3L 

36 

3L 

-L 2 

-3L 4L 
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where the terms within the matrix are the same as the 
consistent geometric stiffness matrix. Thus, [Ro] lacks 
rigid body rotation capabilities. 


CHAPTER 5 


FORCE UNBALANCE 

Closer examination of the traditional static 
formulation of [Kg] indicated that there is a load imbalance 
in the representation, and that pseudo-forces occur to 
maintain equilibrium (Figure 3). 

Recall that [Kg]*{U RBR } = { -P9 , 0 , P0 , 0 } . Using 

Figure 3, and letting the sum of the moments at 0 equal 
zero, yields 

PLSIN2B - P’LCOSB = 0 
P' = P-TAN2B 
= P-TAN0 

= P0 + higher order terms 

Thus, P' represents pseudo-forces required for equilibrium. 

In reference [15], Collar and Simpson acknowledge 
the lack of rigid body rotation capability of [Kg], but 
indicate that it is not a problem, because the energy 
representation is correct. 
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p 


Figure 3 



P' Represents Pseudo-forces Required for 


Equi librium 
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Consider the work/energy relationship from Figure 
3 , without P ' . 

WORK DONE BY P = PL(l-COS2B) 

= 2PL(l-COS2fi)/2 

= 2PL • S IN 2 B 

2 4 

= 2PLB + 0(B ) + higher order terms 

Similarly, using a matrix development 

ENERGY = 1/2{U} T - [K] • {U} 

= PB 2 /2[-2, 0,2,0] • [-L,2,L,2] T 
= 2PLB 2 . 

Therefore, the energy relationship is correct for 

2 

the B terms, but the higher order terms are neglected. For 
large rigid body rotation, this is significant. 

It should be noted that as long as the pre-load P is 
assumed to remain horizontal during rotation, work will be 
done by the force. Thus, true rigid body rotation cannot 
occur. In order for the true strain energy to equal zero, 
the force P must change its orientation as the beam rotates 
(ie. a follower force, as in Figure 4). 

WORK DONE = -L ( P+P • COS2B) ( 1-COS2B ) / 2 + P • SIN2B ( L • SIN2B ) / 2 
= PL[-(l+COS2B) (l-COS2B)+SIN 2 2B] 

= PL/2(-l+COS 2 2B+SIN 2 2fl) 

= PL/2 ( -1+1) 


0 . 
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CHAPTER 6 


RIGOROUS SOLUTION OF FREE/FREE BEAM WITH AXIAL TENSION LOAD 

The author [16] also developed the rigorous free 
vibration solution of a free/free beam with an axial tension 
pre-load. The equation of motion developed agrees with that 
given by Paz [13] uses to develop his dynamic stiffness 
matrix. Solution of the differential equation is similar to 
that given by Shaker [17]. It was also shown that the 
characteristic equation developed indicated the presence of 
three zero frequencies, and the corresponding rigid body 
modes . 
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CHAPTER 7 


DIRECTED FORCE PROBLEM 

Since, as was shown in Section 4, traditional 
formulations did not satisfy equilibrium conditions during 
rigid body rotation, it was determined that investigation of 
the directed force problem is necessary. 

Consider a beam with axial forces which remain 
directed at the opposite end points (Figure 5). This force 
system can be shown to be conservative. Derivation of 
stiffness matrices for this system was examined, utilizing 
Clough' methodology [7], Saunders' methodology [9], and 
Galerkin's criterion. The first two methods are discussed 
in this section. The last method is discussed in Section 9. 

Clough's Methodology for the 
Directed Force Problem 

Consider the beam in the deformed state shown in 
Figure 6. Summing forces in the X direction and setting 
them equal to zero yields 

P cos 0 = N(x) cos W'(x) + V(x) sin W'(x) 
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P cos 0 V(x) sin W’(x) 

N(x) = - (7.1) 

cos W ' ( x ) cos W ' ( x ) 


For small displacements 


N ( x ) = P cos 9 - V ( x ) W' (x) 

Summing forces in the Y direction and setting them equal to 
zero yields 


P sin 0 = N(x) sin W' - V(x) cos W' 


N(x) 


P sin 0 + 


sin W ' 


V ( X ) cos W' 
sin W' 


(7.2) 


For small rotations, Eq.(7.2) becomes 


N(x) = P sin 0 + V(x) 


W' 


W' 


Equations (7.1) and (7.2) can be rewritten 


N(x) cos W' P cos 0 cos W' 


sin W ' 


cos W' sin W' 


V(x) 


(7.3) 
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N(x) sin W P sin 9 sin W’ 

= + V ( x ) (7.4) 

cos W' sin W' cos W' 


Adding equations (7.3) and (7.4) yields 


COS 

w' 

x 

sin 

w' 

= p 

cos 8 

sin 

8 

sin 

T 

w' 

cos 

w 1 

sin W ' 

cos 

W 

_ 



. 


_ 


- 


(7.5) 


or 


N(x) 


P(cos 8 cos W' + sin 8 sin W* 
sin W' cos W' 
cos W' 2 + sin W ' 2 


sin W' cos W' 


N ( x ) = P(cos 8 cos W' + sin 8 sin W') (7.6) 

During rigid body rotation, W' = 0, and equation (7.6) 

becomes 


N ( x ) = P. 

Clough and Penzien [7] develop the equation for the 
geometric stiffness matrix as 
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L 

Kgij = N ( x ) Hi'(x) Hj’(x) dx (7.7) 

0 

The Hermitian interpolating polynomials, and their 
derivatives, are 

H 2 = 1 - 3 ( x/L ) 2 + 2 ( x/L ) 3 H 5 = 3(x/L) 2 -2 ( x/ 1 ) 3 

H 2 ’ = -6 x/L 2 + 6x 2 /L 3 H 5 ' = 6x/L 2 - 6x*/L 3 

7 37 7 37 (7.8) 

H 3 = X -2x / L + x J /L Z H 6 ’ = -x/L + x^/L^ 

H 3 ' = 1 - 4x/L + 3x 2 /L 2 H 6 ' = -2x/L + 3x 2 /L 2 

Substituting equation (7.6) into equation (7.7) yields 
L 

Kgij = P(cos 6 cos W' + sin 0 sinW') Hi'(x) Hj'(x) dx 
0 L 

= P cos 0 cos W' Hi'(x)Hj*(x) dx 
0 L 

+ P sin ©Jsin W Hi' (x)Hj' (x)dx (7.9) 

0 

For the special case of pure rigid body rotation only 
(W'(0) = W'(L) = 0 / f(x)), equation (7.9) reduces to 

Kgij = P(cos 2 0 + sin 2 0) Hi'(x)Hj'(x) dx 

0 
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= P 


Hj_ ' ( x)H j ' (x) dx 


which yields the consistent K g matrix. 

Therefore, K g using this formulation does not 
possess rigid body rotation capability. It may be necessary 
to include separate interpolating functions for the rotation 
instead of the derivatives of the shape functions. 

Consider the appropriateness of the Hermitian 
polynomials for shape functions. 


u 


{Hi H 4 } 


«1 

u 2 


v 


(h 2 h 3 h 5 h 6 } 


u 2 

u 3 

u 5 

U6 


For rigid body translation 

u = ui = u 2 = constant = ft 
a = Hi fl + H 4 Cl 

= (Hi + h 4 ) a 

But, Hi + H 4 = 1, so the equality is satisfied. 

Similarly, v = u 2 = U 5 = constant = 0 
and u 3 = ug = 0 . 



37 


v = H 2 v + H 5 0 
v = (H 2 + H5) 0 


Since H2 + H5 = 1 , the equality is satisfied. 


8 = dv/dx 


{H 2 ’ H 3 ' H 5 ' H 6 '} 


u 2 

^3 

u 5 

u 6 


For rigid body rotation 

0 = u 3 = ug = constant = 0 
u 2 = _u 5 = “L sin 0/2 

8 = H 2 ' (-L sin 0 ) / 2 + H 3 ’ 8 + H 5 ' (L sin 0)/2 +H 6 ' 8 
= (H 5 ' - H 2 ' ) L sin 0 / 2 + (H 3 ’ + Hg') 0 

A <*. 

For small angles, sin 8*0. Therefore, 

8 = [(H 5 ' - H 2 ) L/2 + H 3 ’ + H 6 ’3 0 

0=0. Thus, the equality is satisfied during rigid body 


rotation. 
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Saunders' Methodology for 
the Directed Force Problem 

Saunders considers a beam with a horizontal load as 
shown in Figure 7, and applies various unit displacements 
to develop the stiffness matrices. This can be changed to a 
directed force problem by letting 


P > P = P cos 0 


and 


Vi * v ! = Vi - P sin 0 


as shown in Figure 8. 

S ( x ) = - P cos 0 dy/dx - P sin 0 (7.10) 

= V - P cos 0 y' 

M(x) = V^x + P cos 0 (y^ - y x ) - P x sin 0 - Mi --(7.11) 
= Vix + P cos 0 (yi - y x ) - Mx 

Applying Saunders' methodology using equations (7 . 10 ) and 
(7.11) for the general directed force problem yields 
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M X /EI = 0 ' x (7.12) 

y'x = ® x + r x ( 7.13) 


1 x 

\ 

Rotation due to shear strain 


r x = -6 X /K ' AG (7.14) 

S x = V^- P cos 0 y ' (7.15) 


Therefore , 


r x = V]_ P cos 0 y' 

_ + 

K ' AG K ' AG 


Vi P cos 0(6 X + r x ) 

r x = + 

K ' AG K ' AG 


P cos 0 6 x P cos 0 r x 

r x = + + 

K ' AG K ’ AG K ' AG 


P cos 0 

r x 1 ; 

K AG 


P cos 0 0 X 

+ 

K' AG K’ AG 
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Let R = 1-P cos 0/K'AG 

~ 2 

Vl P cos 0 0 X 

r x + (7.16) 

K ’ AGR K ’ AGR 

Differentiate equation (7.16). 

P cos 0 0 x ' 

r x’ = ; 

K' AGR 

But K ' AGR = (K'AG) (1 - P cos 0 / K ' AG) 

= K'AG - P cos 0. 

Thus , 

P cos 0 0 X ' 

r x' = 

K’AG - P cos 0 

Let P = P cos 0. 

Then , 

r x = P 6 x ' / ( K ' AG - P) 

Differentiate Equation (7.13). 

y x " = + r x ’ 

= 0 X ' + P 0 x ' / ( K ' AG - P) 

= 0 X ' ( 1 + P/ (K'AG - P) ) 


(7.17) 
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Substituting Equation (7.11) into Equation (7.12) yields 


Vjx - Mi + P cos 0 (yi - y x ) 


El 


= ®x' 


Substitute into Equation (7.17). 


v l* - + P(yi-y x )r - -l 

y x " = 1 + P/(K'AG - P 

El L J 


But 


( 1 + P/K'AG - P) = (K'AG - P + P ) / ( K ' AG - P) 

= 1/ ( 1 -P/K'AG) 


Therefore , 


Vix - Mi + P(yi - y x ) 


y x 


(1/(1 - P/K'AG)) 


El 


yx" + p y x / EIR = (Vix - Ml + P yi) / EIR 
Let B 2 = P / EIR 

yx + B 2 y x = B 2 yi - Mi B 2 /P + Vi B 2 x/P (7.18) 

The solution for Equation (7.18) is 


y x = C sin Bx + D cos Bx + yi - Mi/P + Vix/P (7.19) 
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Differentiate Equation (7.19). 


y x ' = CB cos Bx - DB sin Bx + V^/P 


Employing Equations (7.13), (7.14), and (7.15) yields 


®X = Yx' " r x 

= y x ' + S x / K ' AG 

= Yx' + (Vi - P y')/K'AG 
Substitution yields 

0 X = C B cos Bx - D B sin Bx + V^/P 

+ [ Vl - p[c B cos Bx - D B sin Bx + 

K ' AG *- L 

8 X = C B( 1-P/K' AG) cos Bx - D B(1-P/K'AG) sin Bx + 
0 X = C B R cos Bx - D B R sin Bx + Vj/P 

Let y^ = 1 , and apply boundary conditions. 


y ( o ) = yi = l 
0 ( 0 ) = 0 ! = 0 
y ( l ) = y 2 = o 
0 ( l ) = 0 2 = o 


Vi/p] 

Vi/P 

(7.20) 


( See Figure 9 . ) 
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Substitute 9(0) = 0 into Equation (7.20). 

0 = C B R + Vi/P 

C = - Vi/(P B R) 

Substitute 0(L) = 0 into Equation (7.20). 

0 = C B R cos BL - D B R sin BL + Vi/P 
D B R sin BL = C B R cos BL + Vi/P 

C cos BL Vi 

D = + -x 

sin BL P B R sin BL 


Vi cos BL Vi 

D = - x +~ 

P B R sin BL P B R sin BL 


D 


Vj 

PB R 


1 - cos BL 


sin BL 


Similarly, solve equation (7.19) at x = 0 . 

yi = D - Mi + yi 
or 


(7.23) 


D = Mi/P 


(7.25) 
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Solve Equation (7.19) when x = L. 

0 = C sin BL + D cos BL + yi - M^/P + V^L/P 
Yi = Mi/P- ViL/P- C sin BL - D cos BL 
Substituting Equation (7.25) 


yi = D - D cos BL - ViL/P- C sin BL 

7 ! = D ( 1 - cos BL) - ViL/P- C sin BL (7.26) 


Substituting the values for C and D from Equations (7.23) 
and (7.22) yields 


Vl 


Y 1 = ~- 


P B R 


1 - cos BL 


sin BL 


- cos BlJ ViL Vi sin BL 
P PB R 


Vl 


yi = 


(1 - cos BL) sin BL 

L + 


B R sin BL 


B R 


(7.27) 


yi = 


Vl 

2(1 " 

cos BL) - L B R sin BL 

p 


B R sin BL 


Let z 


2(1 - cos BL) - L B R sin BL 
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m* 

Vi Z 

n = » 

P B R sin BL 


Vl - - 

Rll = — = V 1 P B R sin BL / 

yi 

K n = fB R sin BL) / z 

Mi P Vx (1 - cos BL) PB R sin BL 

k 21 = — = * 

Yl (P B R sin BL v i z) 

k 21 = P (1 ‘ cos BL)/z. 

Let 0^ = 1 and apply boundary conditions. (See Figure 10). 


y(0) = yi = 0 

0 ( 0 ) = 0 ! = 1 

Y (L ) = Y2 = 0 

e(L) = 0 2 = o 

Yl = 0 = D - Mi/P (7 . 28) 

Therefore, D = Mi/P (7.29) 


y 2 = 0 - C sin BL + D cos BL - Mi/P + Vix/P 
- C sin BL + D ( 1 - cos BL) = Y L / p 


(7.30) 
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P COS 0 


Figure 10 


Beam with Unit Rotation 0i and 
Other Degrees of Freedom Fixed 
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0 2 = o = C B R cos BL - D B R sin BL + Vj^/P (7.31) 

Therefore, C B R cos BL - D B R sin BL = -V^/P (7.32) 


Solve Equations (7.30) and (7.32) simultaneously. 


- sin BL 

1 - cos BL 


C 


Vi L/P 

B R cos BL 

- B R sin BL 


D 


- Vi L/P 


L V 1 cot BL L Vi Vi 

c = W tp + TP (7.33) 

P P sin BL P B R 


Vl 


BLR sin BL - 1 + cos BL 


P B R 


1 - cos BL 


(7.34) 


D = 


(BLR cos BL - sin BL) 
P B R (cos BL - 1) 


Vl 

D = - tp 

P B R 


BLR cos BL - sin BL 


1 - cos BL 


(7.35) 


Using Equation (7.20) we get 


0! = C B R + Vx/P 


® 1 = — *»" 


Vl 


P B R 


BLR sin BL - 1 + cos BL 


1 - cos BL 


Vl 


Vl 

35 " 

P 


1 - cos BL - B L R sin BL + 1 - cos BL 


1 - cos BL 


6l 



2 cos BL - B L R sin BL 


1 - cos BL 


Vi z 

0 i = -B 

P (1 - cos BL) 


K 12 = Vi/0i 

k 12 = P (1 " cos BL ) / z 


Equations (7.29) and (7.35) yield 


Mi = Vi B L R cos BL - sin BL] 

P P B R 1 - cos BL 


Mi = 


Vl 


BLR cos BL - sin BL 


B R 


1 - cos BL 
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Vl 

BLR cos 

BL - sin BL 


P (1 - cos BL) 

B R 

1 - 

cos BL 


Vi z 


*22 


P sin BL 


z B R 


L cos BL 


*44 


The reaction can be obtained using static equilibrium. 


V 2 = - V 1 (7.36) 

- Mi - M 2 + Vi L + P ji - P y 2 = 0 (7.37) 

or 

M 2 - - M]_ + Vi L + P y] * P 72 (7.38) 


*31 = V 2 /yi = - V]_/yj_ = - P(B R sin BL)/z 
K 31 = - P (B R sin BL)/z = K 13 
and 

M 2 Mi + Vi L + P yi - P Y 2 
*41 = — = 

yi yi 

P(1 - cos BL) P L (B R sin BL) P z 

= + + 

z z z 

= P(-l + cos BL + B L R sin BL + 2 - 2 cos BL 


-BLR sin BL)/z 
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K41 = P (1 - cos BL)/z = K^4 

k 32 = v 2 /e 2 = - Vi/0! 

k 32 = ~ P (1 - cos BL)/z = K 2 3 


Using Equation (7.39) yields K4 2 . 

K 42 = M 2 /0! = (- M X + Vi L + P yi 
P sin BL 

= — — + L cos BL + L ( 1 

z B R 


*42 



sin BL 


B R 


*2 4 


Similarly, 

m 

K33 = P(B R sin BL)/z 
and 

*34 = ~ P(1 " cos BL)/z = K43 


P Y2 ) /®1 
cos BL)] 


or, the final directed stiffness matrix is 
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BR sin BL 1-cos BL - BR sin BL 

sin BL 

1-cos BL Lcos BL cos BL - 1 

BR 


1-cos BL 

sin BL 

L 

B R 


-BR sin BL 
1-cos BL 


cos BL-1 

sin BL 

L 

B R 


BR sin BL cos BL-1 

sin BL 

cos BL-1 — Lcos BL 

BR 


By inspection, K has rigid body translation capability. 
For rigid body rotation (Figure 11), 


U RBR 


L L 

- e, e, - e, e 

2 2 


or 


= [ -L , 2, L, 2] 

K 2j ' U RBR * 0 

K lj ' U RBR =-26 cos BL - B L R B sin BL + 2 0 
= 2 0 (1 - cos BL)) - B L R 0 sin BL 
f 0, No good. 

This matrix does not possess rigid body rotation capability. 

By observation, this [K] is identical with Saunders 
[K] , except P has been replaced with P cos 0, which is the 
component of P in the horizontal direction. It does not 
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possess rigid body rotation capability. It should also be 
noted that the approximation for shear was used to 

determine the stiffness coefficients. Without that 
approximation the stiffness coefficients would be a function 
of the shear. 


ie . 


- Vj P B R sin BL 

Kll = 

(Vi - P sin 0) z 


P (B R sin BL) 


Lim 


z 


= Kll 



Figure 11 Rigid Body Rotation 
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CHAPTER 8 


MATRIX DEVELOPMENT USING GALERKIN CRITERION 

Consider a simply-supported beam with an axial load. 
The differential equation for static displacement is 

E I d 4 W/dx 4 = 0. 


Let W = X 0i(x) Wi 

\ 

Shape functions 


Nodal displacements 


Application of Galerkin's criterion yields 
x 2 

d’W 

0 


x 2 4 

r d 4 w 

El 0j(x) - 7 - dx = 

J dx 4 


Integrate by parts 
Let U = 0 j 


dV = 


d 4 W 


dx 


dx 


d0j 

dU = dx 

dx 


V = 


d 3 W 


dx' 
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3 x 2 
d w 

El 0j (x) 

dx 3 xi 

Integrate by parts. 
Let U = 0j’ 


x 2 3 
f d0 j d W 

EI — 

J dx dx J 


dx = 0 




dv = W 


dU s 0j" dx 


V = 

W" 

3 J 

'2 

d0j 

2 x 2 x : 

d v* 


d w r 

EI 0j(x) — =• 

- 

EI — - 

2 + 

dx 3 


dx 

dx 2 j J 


dx 


X 1 


*1 X 1 


The first two terms are part of the boundary conditions. 
The last term produces the elastic stiffness matrix. 


L 

[K e ] = j EI 0i" 0 j" dx 

0 



H 2 m H2" 

H 3 "H 2 m 

H 3 "H 3 " 

symmetric 

H5 ,, H2 h 

h 5 ,, h 3 " 

h 5 ”h 5 " 

h 6 ,, h 2 m 

h 6 "h 3 " 

h 6 "h 5 " h 6 "h 6 " 


dx 


If the Hermitian polynomials are chosen for the shape 
functions, the resultant matrix is 
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[Ke] = El 


0 

0 12/L 3 

0 6/L 2 

0 0 

0 -12/L 3 

0 6/L 2 


4/L 

0 

-6/L 2 

2/L 


0 

0 

0 


Symmetric 

12/L 3 

-6/L 2 4/L 


This [Ke] has four lero eigenvalues. Note that the K^, 
k 14' k 41' anc * k 44 terms are zero, since we did not develop 
the relationship for the axial terms. Development of these 
terms using a standard Bernoulli formulation yields the 
appropriate terms. If the AE/L terms are developed via 
a classic Bernoulli formulation, the resulting matrix has 
three zero eigenvalues. 


For an axially-loaded beam, the static deformation equation 
becomes 



d 4 w 

d 2 W 

El 

— + P 

dx* 

7? 


Application of Galerkin's criterion yields 
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El 0 


d 4 W 

3 17 


+ P 0 


d 2 w 


dx = 0 


The first term produces [K e ] as shown previously, 
Integrate the second term by parts. 


Let 


[Kg] 


= p J 


u = 

0 


dv = W" 


dU = 

0 

' dx 


V = W' 


L 



L L 


r 

d W 

dW 

r 

dW 

P 0j 


r dx = P 

0 j — 

- P 0 

j’ — dx 

J 

dx' 

i 

dx 

J 

dx 

0 



d 

0 


cond term on the right hand 

side pr 

ic stiffness matrix. 





J 




[Kg] = 

p 

*i' #j' 

dx 









Hi’Hi* 






H 2 , Hi' 


H 2 'H 2 ' 


symmetric 

Hs'Hl’ 


H 3 'H 2 ' 

H 3 'H 3 1 



H 4 ’Hi’ 


H4 , H 2 ' 

H 4 ’H 3 - 

H 4 , H 4 f 


Hs’Hi' 


H 5 'h 2 ' 

KS'H3' 

h 5 , h 4 ' 

H5'H 5 ' 

He'Hl* 


H 6 ’H 2 ’ 

He'Hs’ 

h 6 ’h 4 ' 

H6'H5’ 


dx 
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Choosing the Hermitian Polynomials for the shape functions 
yields 


[Kg] = P 


1/L 

1/L 

6/5L 


0 

1/10 

2L/15 

-1/L 

-1/L 

0 

-1/L 

-6/5L 

-1/10 

0 

1/10 

-L/30 


symmetric 

1/L 

1/L 6/5L 

0 -1/10 2L/15 


Recall that the Ki , j / K 4 , j *Ki , 1 » K i , 4 were all zero 
in the consistent geometric stiffness matrix. When 1/2 
angle of rotation equals one radian, this [Kg] agrees with 
the modified [Kg] developed by Bosela in [2]. Hence, it 
possesses the required three zero eigenvalues. It should 
also be noted that these new terms are not directly 
attributable to the differential equation. 

It should be noted that the consistent [K e ] matrix, 
which has three zero eigenvalues, is utilized along with the 
modified [Kg], which also has three zero eigenvalues, in an 
equation of the form 



a 2 [M] 


= 0 , 


There are only two zero eigenvalues, corresponding 
to rigid body translations. Apparently, the third zero 
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eigenvalue, corresponding to rigid body rotation, is lost, 
due to differences in the rotation eigenvector produced. 

For example, consider the free free beam shown in 
Figure 12. When E= 0 (so that [Ke] is not calculated, and 
only [K] = [Kg] is assembled) and P= 10xl0 6 LBS, the 

following results are obtained (using the computer program 
MODFINITE . FOR : 

Lambda (1) = -0.0010 

Omega (1) = 0.0000 RAD/S 

The associated eigenvector is: 

0.1000000000D+01 
-0. 1000000020 D+ 01 
0. 200000003 9D-01 
-0.9999999973D+00 
0.9999999829D+00 
0.2000000039D-01 

Lambda (2) = 0.0000 

Omega (2) = 0.0000 RAD/S 

The associated eigenvector is: 

0.1000000000D+01 
0. 1166342998 D+ 01 
0.6138287495D-09 
0 . 9999999386D+00 
0.11663430 60D+01 


0. 6138287 603D-09 
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P 


1 

o- 


2 

o 

100 ’ 


P 


A = 48 in 2 

I = 1000 in 4 

m = 0.03525 lb-sec 2 /in 2 


Figure 12 Free/Free Beam with Axial Compression Load 



0.0000 


Lambda (3) 
Omega (3) 


0.0000 RAD/S 


The associated eigenvector is: 


0.1000000000D+01 
-0. 8573806678 D+ 00 
-0.5222345200D-09 
0.1000000052D+01 
-0.8573807200D+00 
-0.5222345360D-09 

Lambda (4) = 

Omega (4) = 


552496.2879 

743.3009 RAD/S 


The associated eigenvector is: 


0.1000000000D+01 
0.6229578457D+00 
0.2524505863D+01 
-0.1000000000D+01 
-0.6229578457D+00 
0. 25245058 63D-01 


Lambda (5) = 1702127.6596 

Omega (5) = 1304.6561 RAD/S 

The associated eigenvector is: 


-0.1266228240D-15 
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0 .1000000000D+01 
-0 . 6000000000D-01 
0.1342393452D-15 
0 .1000000000D+01 
0. 6000000000D-01 


Lambda (6) = 4894312.1933 

Omega (6) = 2212 . 3092 .RAD/S 

The associated eigenvector is: 


0 .1000000000D+01 
0.1337704207D-02 
-0 .1505245050D-01 
-0 . 1000000000D+01 
-0 .1337704207D+02 
-0 . 1505245050D+01 


It can also be shown that the eigenvectors are 
combinations of the rigid body translations, 

[10010 0] T 
[01001 0] T 

T 

and [L -L 2 -L L 2] , which represents rigid body r 
with 1/2 angle or rotation equal to one radian. 

Assuming P=0 and E= 30x10^ psi, the fol 
results (once again using MODFINIT . FOR) are obtained: 


linear 


otation 


1 owing 



Lambda (1) = 0.0000 

Omega (1) = 0.0000 RAD/S 

The associated eigenvector is: 

0.1000000000D+01 

0.0000000000D+00 

0.0000000000D+00 

0.1000000000D+01 

0.0000000000D+00 

O.OOOOOOOOOOD+OO 

Lambda (2) = 0.0000 

Omega (2) = 0.0000 RAD/S 

The associated eigenvector is: 

0.0000000000D+00 
0 . 9334669755D+00 
0. 665330244 6D-03 
0.0000000000D+00 
0.1000000000D+01 


0. 665330244 6D-03 


Lambda (3) = 0.0000 

Omega (3) = 0.0000 RAD/S 

The associated eigenvector is: 
0.0000000000D+00 
0.1000000000D+01 
-0.1977319320D-01 
O.OOOOOOOOOOD+OO 
-0.9773193204D+00 
-0.1977319320D-01 

Lambda (4) = 6127659.5745 

Omega (4) = 2475.4110 RAD/S 

The associated eigenvector is: 
0.0000000000D+00 
0.1000000000D+01 
-0.6000000000D-01 
0.0000000000D+00 
0 . 1000000000D+01 
0 . 6000000000D-01 

Lambda (5) = 49021276.5957 

Omega (5) = 7001.5196 RAD/S 

The associated eigenvector is: 
0.1000000000D+01 
O.OOOOOOOOOOD+OO 


0.0000000000D+00 
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-0 .1000000000D+01 
0.0000000000D+00 
0 . OOOOOOOOOOD + OO 

Lambda (6) = 71489361.7021 

Omega (6) = 8455.1382 RAD/S 

The associated eigenvector is 
0. OOOOOOOOOOD + OO 
0.1000000000D+01 
-0.1200000000D+00 
0 .0000000000D+00 
-0.1000000000D+01 
-012000000000 D+ 00 

Once again, the first three eigenvectors can be 
shown to be linear combinations of rigid body translational 
modes . 

[100 10 0] T 
[010 01 0] T 

and the rigid body rotation vector 
[ 0 -L 2 0 L 2] 

Considering the same beam, with P=10X10 6 lbs and 
E=30X10 6 , the following results are obtained from 


MODFINIT . FOR . 



|[Ke] + P[Kg] - Q^[M]| = 0 


Lambda (1) 
Omega (1) 


-0.0688 
0.0000 RAD/S 


The associated eigenvector is: 
0.1000000000D+01 
-0. 10000000 16D+01 
-0. 157767551 4D -15 
0.1000000000D+01 
-0.1000000016D+01 
-0.3281352645D-15 


Lambda ( 2 ) 
Omega ( 2 ) 


0.0688 

0.2622 RAD/S 


The associated eigenvector is: 
0.1000000000D+01 
-0.999999984 1D+00 
-0.1577674314D-15 
0 . 1000000000D+01 
-0 . 9999999841D+00 
-0.3281350580D-15 


Lambda (3) 
Omega (3) 


403295.7555 

635.0557 RAD/S 


The associated eigenvector is: 


0 . 1000000000D+01 



-0. 1066270395 D+ 01 


0.1239429315D-02 
0.100018382 9D+ 01 
-0. 9339134336D+00 
0.1239429315D-02 

Lambda (4) = 7829787.2340 

Omega (4) = 2798.1757 RAD/ S 

The associated eigenvector is: 

0 .1000000000D+01 
-0.3571250023D+00 
-0.3857249986D-01 
0.1000000000D+01 
-0.3571250023D+00 
0. 385724998 6D -01 

Lambda (5) = 76419084.4813 

Omega (5) = 8741.8010 RAD/S 


The associated eigenvector is: 
0.1000000000D+01 
-0. 7447662297 D+ 01 
0.7700947165D+00 
0 .1017885331D+01 
0 . 5429776966D+01 


0.7700947165D+00 
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Lambda (6) = ******************** 

Omega (6) = 248742850249.1149 RAD/S 

The associated eigenvector is: 

0.1000000000D+01 

-0.1000000000D+01 

0.2272727273D-01 

-0.1272727273D+01 

0.1272727273D+01 

0.2272727273D-01 

If the [K e ] matrix generated by the Galerkin method 
(which has 4 zero eigenvalues) is used along with modified 
[Kg] from Galerkin (which has three zero eigenvalues), the 
combined stiffness yields three zero eigenvalues. 

The results are obtained from BOFINITE . FOR : 

| [ [Ke] + P[Kg] - Q 2 [M]]| = 0 

yields 

Lambda (1) = -0.0010 

Omega (1) = 0.0000 RAD/S 

The associated eigenvector is: 

0.1000000000D+01 
-0. 9999995972D+00 
0.1999999995D-01 


-0 . 9999999870D+00 



0.1000000396D+01 


0.1999999995D-01 

Lambda (2) = 0.0000 

Omega (2) = 0.0000 RAD/S 


The associated eigenvector is: 
0.1000000000D+01 
0.1686129333D+02 
-0.2021558453D-06 
0 . 100000039 6D+01 
0.1999999995D-01 


Lambda ( 2 ) 
Omega (2) 


0.0000 

0.0000 RAD/S 


The associated eigenvector is: 
0 . 1000000000D+01 
0.1686129333D+02 
-0.2021558453D-06 
0.1000020216D+01 
0. 1686127311 D+ 02 
-0.2021558449D-06 


Lambda (3) 
Omega (3) 


0.0000 

0.0000 RAD/S 


The associated eigenvector is: 


0.1000000000D+01 



-0.5930809082D-01 


0.5153549492D-09 
0.9999999485D+00 
-0.5930803928D-01 
0.5153549432D-09 

Lambda (4) = 

Omega ( 4 ) = 

The associated eigenvector is: 

0.1000000000D+01 
0. 9739809665 D+ 00 
-0.1747771626D-01 
-0.1000000000D+01 
-0. 9789809665D+00 
-0.1747771626D-01 

Lambda (5) = 7829787.2340 

Omega (5) = 2798.1757 RAD/S 

The associated eigenvector is: 

0. 280148715 6D-16 
0.1000000000D+01 
-0.60 0000000 0D -01 
-0.2859899241D-16 
0.1000000000D+01 


673695.6482 
820.7896 RAD/S 


0 . 6000000000D-01 
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Lambda (6) = 76262474.5352 

Omega (6) = 8732.8389 RAD/S 

The associated eigenvector is: 


0.1000000000D+01 
0.2230210189D+03 
-0.2666252228D+02 
-0.1000000000D+01 
-0. 223021018 9D+03 
-0.2666252228D+02 


Next consider the equation of motion 

d 4 W d 2 W d 2 W 

El — t + N — =■ + m — ~ = 0 
dx* dx^ dr 

The first two terms have already been examined. 

2 2 

consider only the term m d W/dt . 

Let W = I di(x) Wi 
and 


W(x,t) = 0(x) sin Qt 
(d/dt) W(x,t) = 0 ( x ) Q cos Qt 
(d 2 /dt 2 ) W(x,t) = - 0 ( x ) Q 2 sin Qt 
= - Q 2 W(x,t) 


Now 


Apply Galerkin's criterion 
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x 2 L 

| m 0i d 2 W/dt 2 = 6 A | 0i 0j dx 
xi 0 


[H] 



Hl«i 


symmetric 


H 2 H! H 2 H 2 


H 3 H 1 

H 3 H 2 

H 3 H 3 



H 4 H 1 

H 4 H 2 

H 4 H 3 

H 4 H 4 


H 5 H 1 

H 5 H 2 

H 5 H 3 

H 5 H 4 

H 5 H 5 

H 6 Hi 

H 6 H 2 

H6«3 

H 6 H 4 

H6«5 



Once again, selecting the Hermitian polynomials for the 
shape functions yields 


140 

147 156 symmetric 

SAL 21L 22L 4L 2 

[M] = 

420 70 63 14L 140 

63 54 13L 147 156 

-14L -13L -3L 2 -21L -22L 

Note that the M 12 , M 13 , M^, Mjg, etc. are zero in 
the consistent mass matrix, and not directly attributable to 
the differential equation used. 

In order for finite element dynamic analysis 
alorithims to provide solutions, the mass matrix must always 
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be positive definite. To test whether the mass matrix would 
be positive definite, consider the following: 

Let m = 0.03525 LB Sec 2 /IN 2 
L = 10 IN 

The mass matrix becomes 


0.1175 

0.1234 

0.1762 

0.05875 

0.05287 

-0.1175 

0.1234 

0.1309 

0.1846 

0.05287 

0.04532 

-0.1091 

0.1762 

0.1846 

0.3357 

0.1175 

0.1091 

-0.2518 

0.05875 

0.05287 

0.1175 

0.1175 

0.1234 

-0.1762 

0.05287 

0.04532 

0.1091 

0.1234 

0.1309 

-0.1846 

0.1175 

-0.1091 

-0.2518 

-0.1762 

-0.1846 

0.3357 

A positive definite 

ma t r ix has all 

posi tive 


eigenvalues. Solving the algebraic eigenvalue problem, 

| CM] - [I] J =0 yields 

Eig . ! = -3.50 x 10 * 5 
Eig . 2 = 1.41 x 10” 2 

Eig . 3 = 9.02 x 10 * 

Eig . 4 = -4.61 x 10 ‘ 5 
Eig . 5 = 2.14 x 10 -1 

Eig. g = 3.83 x 10 -2 

The negative eigenvalues indicate that this particular mass 
matrix is not positive definite. 



CHAPTER 9 


SAMPLE PROBLEMS 

Fertis and Lee [19] developed the equations of 
motion and obtained rigorous solutions for beams with 
various loading and end conditions. The Galerkin method can 
be used to generate stiffness and mass matrices for a finite 
element application, and the results compared with Fertis 
and Lee's rigorous solution. 

The following beams were considered: 

a. Axially-loaded beam on simple supports. 

b. Axially-loaded beam with vertical spring supports. 

c. Axial 1 y- 1 oaded beam with horizontal and vertical 
support springs. 

d. Bow-string 

It should be noted that Fertis and Lee's analysis 
indicate regions of dynamic "flutter" instability, which has 
not previously been identified. Kounadis [20] has 
identified similar areas in the stability analysis of beams 
with follower forces. 
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Case a. Axial 1 y- 1 oaded Beam on Simple Supports 

Fertis and Lee [19] develop the general equation of 
motion in the form 


dx‘ 


EI(x) 


d*y 

dx 2 


d 2 y 


d 2 y 


+ 6 A( x) — 5 - P — =■ 
dr dx^ 


d 

6 — 
d s 


,3 

d y 

x 2 
dt ax 


(9.1) 


If a beam with a constant cross-section is considered, EQ 
(9.1) reduces to 


d 4 y d 2 y d 2 y d 4 y 

El j + 6 A x - P =■ - 61 = 5 = 0 (9.2) 

dx 4 dt^ dx^ dt ax^ 


It can easily be seen that the first three terms of 
EQ (9.2) yield K E , H, and K g matrices previously developed 
using Galerkin criterion, (with the negative sign on the 
third term indicative of an axial compression load). 


Consider the fourth term. 


6 2 1 


d 4 y 


K'G dt 2 dx 2 



Let fl = I 0i(x) W A 
and 
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W(x,t)=0(x) Sin Qt 
d 


dt 

,2 


dt‘ 


W = 0 ( x ) Q Cos Qt 


W s -0(x)Q Sin Qt 


dt‘ 


W = -Q^ W ( x , t ) 


— = — =■ W = -Q* — ,W ( x , t ) 
dx at^ dx^ 


Coinparison with term #3 indicates that the following 
higher order matrix is produced: 


- Q' 


1/L 




1/L 

6/5L 


SYMMETRIC 

0 

1/10 

2L/15 


-1/L 

-1/L 

0 

1/L 

-1/L 

-6/5L 

-1/10 

1/L 6/5L 

0 

1/10 

-L/30 

0 -1/10 2L/15 

jr Kil 

, and 

Ki 4 terms 

could be set equal to 


(corresponding to more traditional development) since they 
are not explicitly developed using the differential 
equation. 
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Case b. Axial 1 y- 1 oaded Beam with Vertical Spring Supports 

Fertis and Lee [19] have utilized Hamilton's principle 
and the dynamic equilibrium method to formulate the 
characteristic equations for a beam with an axial 
compression load and vertical spring supports. They have 
solved for the frequencies of vibration for various 
parameters utilizing a bisection method, and tabulated the 
results . 

Using Galerkin's criterion, stiffness and mass 
matrices can be generated, formulated into a finite element 
algorithm, and the finite element solution compared with 
Fertis and Lee's rigorous solution. 

Fertis and Lee's analysis indicates that regions of 
dynamic ("flutter") instabilities occur for an axially- 
loaded beam with spring supports. The modified finite 
element approach provides correlation with the trial and 
error procedure. 

The differential equation developed by Fertis and Lee 

[19] is 


El 


4*4 
d y El 6 d y 

* 

5 I 

K’G d 4 y d 4 y 

£ 

dx* K G dx at^ 

k'g 

dx 2 dt 2 

4 

dt 4 

( 1 ) ( 2 ) 


(3) 

(4) 


+m 


d 2 P d 2 y 
dt 2 dx 2 

(5) ( 6 ) 

(9.3) 


04 



with I 
Term (1 

[K e ] 

Terms ( 

[ r rot3 

Term (4 
[Mi] = 


= 1/ ( 1+P/K ' AG ) 
) yields 


12/L 3 





symmetric 

6/L 

4/L 


-12/L 3 

-6/L 2 

12/L 3 

6/L 2 

2/L 

-6/L 2 


2) and (3) yield 


6EI * 

+ 61* 

K’G 


6/5L 

1/10 

2L/15 

symmetric 

-6/5L 

-1/10 

12/L 3 

1/10 

-L/30 

-6/L 2 


4/L 


) yields 


6 2 I* 


420 K' G 


156 


22L 

4L 2 


54 

13L 

156 

-13L 

-3L2 

-22L 



Term (5) yields 
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SAL 

[M 0 ] = 

420 


156 


22L 

4L 2 

54 

13L 

-13L 

-3L 2 


symmetric 

156 

-22L 4L 2 


Term (6) yields 



6/5L 

1/10 

2L/15 

symmetric 

[Kg] = P 

-6/5L 

-1/10 

6/5L 


1/10 

-L/30 

-1/10 2L/15 


In matrix form, EQ (9.3) yields 


|[K E ] + P[K G ]1 - Q 2 j [M 0 ] - [K R0T ]| + Q 4 [Mi] = 0 
L J L J (9 .4) 

The above matrices were included in a finite element 

4 

dynamic analysis program (NLFIN.FOR), neglecting the W 
term. Tables 2-6 compare the finite element output with 
Fertis and Lee's rigorous solution of the dynamic 
instability loads. 
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AXIAL 

COMP 

LOAD 

VK/M 

RIGOROUS 
1ST FREQ 
RAD/ SEC 

NLFIN3 

1ST 

FREQ. 

% 

DIFF. 

RIGOROUS 
2ND FREQ 

NLFIN3 
2ND FREQ 

% 

DIFF 

0 

0.688 

0.689 

0.6881 

0.01 

1.189 

1.1930 

0.34 

100 

0.688 

0.689 

0.6881 

0.01 

1.063 

1.0670 

0.38 

200 

0.688 

0.689 

0.6881 

0.01 

0.921 

0.9241 

0.34 

300 

0.688 

0.689 

0.6881 

0.01 

0.752 

0.7545 

0.33 

355 

0.688* 

0.689* 

0.6881 

0.0 

0.689* 

0.6881 

0.13 


♦Flutter (dymanic instability) occurs. 


















Table 4 Rigorous Solution Versus NLFIN3 . FOR 


A = 
E i 
I s 


AXIAL 

COMP 

LOAD 


K=100 LB/IN 
■ 48 IN 2 
' 30 x 10 6 
* 256 IN 4 


RIGOROUS 
VK/M 1ST FREQ 
RAD/ SEC 


2.176 2.177 
2.176 2.177 
2.176 2.177 
2.176 2.177 
2.176 2.177 
2.176 2.177 
2.176 2.177 
2.176 2.177* 


K=100 LB/IN 


m = 0.4224 LB-SEC 2 /IN 2 


L = 100 IN 
K ' = 0.186 


NLFIN3 

1ST 

FREQ 


2.176 

2.176 

2.176 

2.176 

2.176 

2.176 

2.176 

2.176 


RIGOROUS 
2ND FREQ 

NLFIN3 
2ND FREQ 

% 

DIFF 

3.758 

3.7725 

0.39 

3.565 

3.5789 

0.39 

3.362 

3.3742 

0.36 

3.144 

3.1563 

0.37 

2.911 

2.2922 

0.38 

2.658 

2.6676 

0.36 

2.377 

2.3859 

0.37 

2.177* 

2.1671 

-0.45 


* Flutter (dynamic instability) occurs 



















86 


Table 5 Rigorous Solution Versus NLFIN3 .FOR 
P P 



K=500 LB/IN K=500 LB/IN 

A = 48 IN 2 I * 256 IN 4 m = 0.4224 LB-SEC 2 /IN 2 

E = 30 x 10 6 K'= 0.186 L = 100 IN 

AXIAL RIGOROUS NLFIN3 % RIGOROUS NLFIN3 % 

COMP VK/M 1ST FREQ 1ST FREQ DIFF 2ND FREQ 2ND FREQ DIFF 


LOAD RAD/ SEC 

0 4.866 4.877 4.863 

1000 4.866 4.877 4.863 

2000 4.866 4.877 4.863 

3000 4.866 4.877 4.863 

4000 4.866 4.877 4.863 

5000 4.866 4.877 4.863 

6000 4.866 4.877 4.863 

7000 4.866 4.877 4.863 

8000 4.866 4.877 4.863 

9000 4.866 4.877 4.863 

10000 4.866 4.877 4.863 

11000 4.866 4.877 4.863 

12000 4.866 4.877 4.863 

13000 4.866 4.877 4.863 

14000 4.866 4.877 4.863 

15000 4.866 4.877 4.863 

16000 4.866 4.877 4.863 

16550 4.866 4.877* 4.863 


0.06 

8.402 

8.4350 

0.39 

0.06 

8.232 

8.2646 

0.40 

0.06 

8.057 

8.0907 

0.42 

0.06 

7.882 

7.9128 

0.39 

0.06 

7.701 

7.7309 

0.39 

0.06 

7.510 

7.5446 

0.46 

RKTS 

7.325 

7.3536 

0.39 


7.130 

7.1575 

0.39 


6.929 

6.9559 

0.39 

0.06 

6.722 

6.7482 

0.39 

0.06 

6.509 

6.5339 

0.38 

0.06 

6.288 

6.3124 

0.39 

0.06 

6.059 

6.0828 

0.39 

0.06 

5.822 

5.8422 

0.38 

0.06 

5.315 

5.3350 

0.38 

0.06 

5.315 

5.3350 

0.38 

0.06 

5.024 

5.0612 

0.74 

0.06 

4.877* 

4.9042 

0.56 
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Table 6 Rigorous Solution Versus NLFIN3 . FOR 


K=1000 LB/IN 


48 IN' 


K=1000 LB/IN 

m = 0.4224 LB-SEC 2 /IN 2 


= 30 x 10 


100 IN 


= 256 IN 


K ' = 0.186 


AXIAL 

COMP 

LOAD 

f K/M 

1ST 
FREQ 
RAD/ SEC 

1ST 

FREQ 

DIFF 

2ND FREQ 

2ND FREQ 

DIFF 

0 

6.881 

6.876 

6.874 

0.10 

11.882 

11 . 9282 

0 . 39 

5000 

6.881 

6.876 

6.874 

0.10 

11.272 

11.3162 

0 .39 

10000 

6.881 

6.876 

6.874 

0.10 

10.628 

11.6692 

0.39 

15000 

6.881 

6.876 

6.874 

0.10 

9.941 

9.9802 

0.39 

20000 

6.881 

6.876 

6.874 

0.10 

9.294 

9.2400 

-0.58 

25000 

6.881 

6.876 

6.874 

0.10 

8.402 

8.4350 

0.39 

30000 

6.881 

6.876 

6.874 

0.10 

7.516 

7.5446 

0.38 

31000 

6.881 

6.876 

6.874 

0.10 

7.325 

7.3536 

0.39 

33501 

6.881* 

6.876* 

6.874 

0.10 

6.876* 

6.8526 

-0 . 34 


* Flutter (dynamic instability) occurs. 
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NLFIN3 has excellent correlation with rigorous 
solution (See Tables 3-6). However, the rigorous solution 
has the 2nd frequency varying with axial load, even when 
K=0. Thus, it does not model rigid body rotation. NLFIN3 
was developed from the differential equation in the rigorous 
solution, using Galerkin's criterion. Thus, it also has the 
2nd frequency varying with the axial load, and does not have 
rigid body rotation capability. Table 7 compares the 
critical load obtained using the finite element method 
versus the rigorous solution for varying spring stiffness. 
Tables 8-10 compare the frequencies for constant spring 
stiffness but varying lengths. 

The critical loads correspond to a coalescing of the 
first and second flexural ei genf r equenc i es (dynamic 
instability) and were located by varying the load. 
Correlation between the finite element and rigorous solution 
is evident. The difference is affected by the judgment as 
to when the first two frequencies have sufficiently 

A 

coalesced, the excluding of the Q contribution in the 
finite element algorithm, the number of elements used, and 
errors due to the iterative nature of the eigensolver 
routine . 

It should also be noted that in general, as the 
axial compression load is increased, the second eigenvalue 
decreases. The first eigenvalue does not change 
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Table 7 Critical Load Rigorous Solution Versus NLFIN3 
p o £ P 

A = 48 IN 2 I = 256 IN 4 m = 0.4224 LB-SEC 2 /IN 2 


E = 30 : 

it 10 6 K' = 0.186 

L = 100 IN 


SPRING 

CRITICAL LOAD (LB) 


CONSTANT 

RIGOROUS 

NFLIN3 

DIFF 

10 

3.35 X 10 2 

3.335 X 10 2 

-0.4 

100 

3.35 X 10 3 

3.335 X 10 3 

-0.4 

500 

1.655 X 10 4 

1.6675 X 10 4 

0.8 

1000 

3.100 X 10 4 * 

3.335 X 10 4 

7.6 

2000 

6.600 X 10 4 

6.670 X 10 4 

1.1 

3000 

9.550 X 10 4 

1.0025 X 10 5 

5.0 

5000 

1.605 X 10 5 

1.6708 X 10 5 

4.1 

10000 

3.345 X 10 5 

3.3722 X 10 5 

0.8 

20000 

6.75 X 10 5 

6.8141 X 10 5 

0.9 

60000 

2.105 X 10 6 

2.1523 X 10 6 

2.2 

80000 

2.98 X 10 6 

2.9649 X 10 6 

-0.5 

100000 

3.899 X 10 6 

3.8548 X 10 6 

-1.1 

120000 

4.858 X 10 6 

4.8700 X 10 6 

0.2 

140000 

6.250 X 10 6 

6.1600 X 10 6 

-1.4 

147000 

7.000 X 10 6 

6.7850 X 10 6 

-3.1 

i 

150000 

7.450 X 10 6 

7.1400 X 10 6 

CM 

• 

t 

200000 

7.450 X 10 6 

7.6400 X 10 6 

2.6 

300000 

7.450 X 10 6 

7.6400 X 10 6 

2.6 

400000 

7.450 X 10 6 

7.6400 X 10 6 

2.6 

500000 

7.450 X 10 6 

7.6400 X 10 6 

2.6 
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Table 8 Natural Frequency Versus Length 
P o jp- P 

A = 48 IN 2 I = 256 IN 4 m = 0.4224 LB-SEC 2 /IN 2 

E = 30 x 10 6 K'= 0.666 P = 1,000,000 lb 

Spring constant K = 100,000 LB/IN 


LENGTH RIGOROUS 


NLFIN3 



123.18 123.188 

121.140 121.208 

117.439 117.505 

107.961 108.040 

105.235 105.320 

102.683 102.772 

100.284 100.377 

98.022 98.119 

95.880 95.982 

86.602 86.726 

78.990 79.138 

72.424 72.596 

61.110 61.328 

56.009 56.247 

34.069 34.340 

20.903 21.147 

18.182 18.420 

6.329 6.632 4.8 

3.430 3.892 13.5 

0 1.742 


RIGOROUS 
2ND FREQ 

NLFIN3 
2ND FREQ 

123.18 

128.860 

125.57 

130.254 

127.653 

132.193 

130.602 

134.802 

130.671 

134.148 

130.470 

133.347 

130.042 

132.445 

129.461 

132.437 

128.757 

131.449 

124.222 

125.956 

119.115 

120.344 

114.079 

115.032 

104.831 

105.574 

100.641 

101.376 

82.686 

83.882 

70.114 

72.098 

67.070 

69.288 

52.364 

55.738 

49.577 

53.140 

48.211 

51.860 
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Table 9 Natural Frequency Versus Length 
P 1 o £- P 

A = 48 IN 2 I = 256 IN 4 m = 0.4224 LB-SEC 2 /IN 2 

E = 30 x 10 6 K ' - 0.666 P = 1,000,000 lb 


Spring constant K = 50,000 LB/IN 


LENGTH 

RIGOROUS 
1ST FREQ 

NLFIN3 
1ST FREQ 

% 

DIFF 

RIGOROUS 
2ND FREQ 

NLFIN3 
2ND FREQ 

% 

DIFF 

60.5 

62.023 

61.784 

0.4 

62.023 

63.316 

Bj 

62 

60.977 

60.972 

0.0 

63.197 

63.808 

l 

64 

59.889 

59.930 

0.1 

63.969 

64.546 

0 . 9 

66 

58.885 

58.928 

0.1 

64.585 

65.132 

0.8 

68 

57.919 

57.964 

0.1 

65.070 

65.590 

0.8 

70 

56.987 

57.034 

0.1 

65.445 

65.939 

0.8 

72 

56.087 

56.136 

0.1 

65.726 

66.197 

1 

74 

55.217 

55.268 

0.1 

65.928 

66.376 


76 

54.374 

54.427 

0.1 

66.061 

66.490 

0.6 

78 

53.555 

53.610 

0.1 

66.137 

66.546 

0.6 

80 

52.760 

52.817 

0.1 

66.161 

66.554 

0.6 

82 

51.986 

52.046 

0.1 

66.143 

66.520 

0.6 

84 

51.232 

51.294 

0.1 

66.086 

66.449 > 

0.5 

88 

49.777 

49.843 

0.1 

65.880 

66.218 

0.5 

90 

49.074 

49.143 

0.1 

65.737 

66.064 

0.5 

100 

45.756 

45.838 

0.2 

64.752 

65.082 

0.5 

140 

34.326 

34.462 

0.4 

59.105 

59.418 

0.5 

180 

24.118 

24.296 

0.7 

53.160 

53.651 

0.9 

240 

10.687 

10.910 

2.1 

44.393 

45.492 

2.5 

275 

0 

1.741 

— — 

39.143 

40.768 

4.2 
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Table 10 Natural Frequency Versus Length 
P o o o- P 


A = 48 IN" I = 256 IN’* m = 0.4224 LB-SE 
E = 30 x 10 6 K'= 0.666 P = 1,000,000 lb 
Spring constant K = 25,000 LB/IN 


rr - - 

m = 0.4224 LB-SEC 2 /IN 2 


LENGTH 

RIGOROUS 
1ST FREQ 

NLFIN3 
1ST FREQ 

115.45 

30.529 

30.4991 

120 

29.676 

29.7156 

125 

28.835 

28.8787 

130 

28.015 

28.0628 

135 

27.212 

27.2636 

140 

26.421 

26.4774 

145 

25.640 

25.7010 

150 

24.886 

24.9315 

155 

24.095 

24.1661 

160 

23.326 

23.4026 

165 

22.557 

22.6390 

170 

21.786 

21.8734 

180 

20.231 

20.3305 

190 

18.652 

18.7644 

200 

17.043 

17.1689 

220 

13.718 

13.8721 

i 

230 

11.988 

12.1592 

240 

10.192 

10.3844 

260 

6.164 

6.4413 

275 

0 

1.7385 


RIGOROUS 
2ND FREQ 

NLFIN3 
2ND FREQ 

30.529 

30.6904 

31.269 

31.3525 

31.827 

31.9097 

32.243 

32.3257 

32.547 

32.6295 

32.759 

32.8430 

32.897 

32.9829 

32.974 

33.0624 

33.000 

33.0919 

32.983 

33.0797 

32.931 

33.0326 

32.848 

32.9560 

32.608 

32.7313 

32.291 

32.4328 

31.915 

32.0793 

31.032 

31.2548 

30.539 

30.7994 

30.019 

30.3215 

28.610 

29.3091 

28.211 

28.5057 
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considerably. When the first two frequencies coalesce, 
dynamic (flutter) instability occurs. 

NFLIN3 has excellent correlation with the rigorous 
solution. However, the rigorous solution has the 2nd 
frequency varying with axial load, even when K=0, Thus, it 
does not model rigid body rotation. NLFIN3 was developed 
from the differential equation in the rigorous solution, 
using Galerkin's criterion. Thus, it also has the 2nd 
frequency varying with the axial load, and does not have 
rigid body rotation capability. 

Case c. Axial 1 y- 1 oaded Beam with Horizontal and 
Vertical Support Springs 

This beam was analyzed using a two-element 
formulation, and the results tabulated in Table 11. There 
is excellent correlation between the critical load obtained 
using finite element model and the rigorous solution. 

A four-element model was also used, and the results 
tabulated in Table 12. Once again, there is excellent 
correlation between the frequencies of vibration obtained 
using the finite element model with the rigorous solution. 

It should be noted that for this particular problem, 
flutter does not occur. Instead, instability occurs when 
the natural frequency drops to zero. 


9 


Table 11 Vibration of an Axially-loaded Beam with 
Horizontal and Vertical Springs 



0 

K 

1 

K 



> 

O 

l Kv 



r 

> 

*h 

= 10 x 10 6 LB/IN 

I = 256 IN 4 

K v 

= 400 X 10 3 LB/IN 

E = 30 X 10 6 PSI 

L = 

100 IN 

m = 0.4224 LB-SEC 2 /IN 2 

A = 

48 IN 2 



LOAD 

(LB) 

1ST FREQUENCY 
USING NLFIN5 
(RAD/S) 

0 

93.38 

7,600,000 

9.27 

7,610,000 

7.92 

7,630,000 

4.01 

7,640,000 

2.71 i 

7.635.000 

7.560.000 

2.09 


1ST FREQUENCY 
RIGOROUS SOL 
(RAD/S) 


7,637,000-7,560,000 

7,560,000 


% DIFF 


1.0 % 






Table 12 


Vibration of Axially-loaded Beam 
with Horizontal and Vertical Springs 


Using 4 elements 



X G = 0.0 
P = 0.0 


Freq . 

Rig (RAD/SEC) 

Hodfin 3 

1 

97.8920 

98.1914 

2 

221 . 4823 

224.4996 

3 

385.4009 

419.3890 


X D = 0.225 IN 
P = 2,250,000 LB 


Freq 

Rig (RAD/SEC) 

Hodfin 3 

1 

82.1998 

88.8739 

2 

214.0995 

212.0427 



3 


387.5971 


389.7464 
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Case d._ Bowstring 

Finally, consider a beam which has a spring 
connected between the end nodes, such that the force in the 
spring is always directed between the end nodes. (See Figure 
13 ) . 

The differential equation for the bowstring developed by 
Fertis, is 


#1 #2 #3 #4 #5 #6 



( 9 . 5 ) 


Where 

I* = 1/ (1+P/K' AG) 


K' = Shear Coefficient (shape factor) 
c = Spring Constant 
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Consider term #1 


d 4 y 




El* 

dx 4 



As done previously, applying 

the Galerkin criterion 


12/L 3 

6/L 2 

-12/L 3 

6/L 2 


[Ke ] =EI * 

6/L 2 

4/L 

-6/L 2 

2/L 



-12/L 3 - 

6/L 2 

12/L 3 

-6/L 2 



6/L 2 

2/L 

-6/L 2 

4/L 


Now consider 

the shear 

and rotatory inertia effects 

Consider the 

term 


d 4 y 






, 2,. 2 
dx dt 



Let W(x, 

,t) = Z0(x) 

sin 

Qt 




,4 

d w 

. 2 ji.2" 

dx dt 

-Q 2 

d 2 W 

0(x) 7 

dx^ 
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Apply Galerkin's criterion 


x 2 2 
r d^W 




0-j r dx = 0 

dx^ 


Let u=0j dv = d W 

17 


X 1 


dx 


d0 j 

du = dx 

dx 


v = 


dW 

dx 


x2 


L L 


d 2 W dW 

0j — ^ dx = 0 j — 
dx^ dx 


xl 


dW 

0 j ' — dx 
dx 


x = 0 0 


which yields 


CKRotl 


L 

1 I 

0 i 0 j dx 

0 


Hence, recall term #2. 


El 


+ I 


K G 


d 4 y 


.*2 . 2 
dt dx 


As done previously, applying the Galerkin criterion yields 



100 


[ r rot] = 6 


El 


K'G 


+ I 


6 

5L 

1 

10 

-6 

5L 

1 

10 


1 

10 

2L 

15 

-1 

10 

-L 

30 


-6 

5L 

-1 

10 

6 

5L 

-1 

10 


Where 6 = Mass Density I = 

E = Young's Modulus P = 

I = Moment of Inertia 

K'= Shear Coefficient (Shape Factor) 

G = Modulus of Rigidity 


Consider term #3, 


m 


d 2 y 


As done previously, applying the Galerkin crit 



156 


SYMMETRIC 

SAL 

22L 

4L 2 


420 

54 

13L 

156 


-13L 

-3L 2 

-22L 4L 2 


1 

10 

-L 

30 

-1 

10 

2L 

15 

1/(1+ P/K'AG) 
Axial Load 


erion yields 



Consider term #4, 


d 2 y 

dx 2 


As done previously, applying the Galerkin criterion yields 


[Kg] = P 


6/5L 

SYMMETRIC 


1/10 

2L/15 


-6/5L 

-1/10 6/5L 


1/10 

-L/30 -1/10 

2L/15 


Consider term #6, 


c 2 * ,4 
6 1 d y 


K G dt 


W(x,t) = I0(x)sin Qt 


d 4 W 


dt 


= Q 0 ( x ) Sin Qt 


d 4 W 


dt 


4 

= Q W 


Application of Galerkin's Criterion yields 
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6 1 L 

[Mx] = r~ 

( 420 )K G 


156 


22L 

4L 2 

54 

13L 

-13L 

-3L 2 


Finally, consider the term 


L 

c 

2 


dy 

dx 



dx 


0 


Integrate by parts. 


2 

Let U = y' 


dv 


du = 2 y' y" dx 


v 


L 1 

J ii 

r 2 3 


j y y dx = y ' 

— 

J 

0 ( 

0 


L L 

[ 2 » 3 

3 y' y dx = y' 

0 A 


SYMMETRIC 

156 

-22L 4L 2 


y" dx 
y' 

2 

y" dx 
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Therefore 



Let y = I Yi 

Apply Galerkin’s Criterion 


L L 



0 


3 

Evaluate the term Hi' ^ 
H 2 = 1 - 3(X/L) 2 + 2(X/L) 3 


6x 6x 2 
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= [1 - 4 + 3] 3 -1 
= - 1 . 

H 5 = 3( x/L) 2 - 2 ( x/L ) 3 


6x 6x 2 




105 
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Which yields 

0 0 0 0 

0 -L 2 / 12 0 L 2 /12 

0 0 0 0 

0 L 2 /12 0 -L 2 /12 

By observation, Kg ow possesses both rigid body translation 
and rotation capability. 

For the pre-loaded beam in space, one must consider 
the bow string without the horizontally applied P force, and 
replace the spring force with a constant P force in the 
direction of the spring (Figure 14). 

From Galerkin, the matrix equation included the terms 
C 

P[K g ] - - [K bow ] = 0 
6 

or 

C[Kbow^ = 6 P [Kg] 

One can approximate the condition in Figure 14 by replacing C 
with 6P/L . 
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Thus, one obtains 


E K bow3 


PL 

12 


0 0 

0 -1 

0 0 


0 

0 

0 


0 

1 

0 


0 


1 


0 


-1 


Consider the performance of Kb ow in a buckling 

problem. 

Case 1 Simply-supported Beam (Figure 15). 

This problem is traditionally solved using [Kg] and 
[Kg]. However, due to the nodal restraints in the vertical 
direction, the P forces remain directed at the opposite 
nodes. Thus, [Kg and [K]-, ow ] should provide a similar 
solution. The well-known rigorous solution is 

P = k 2 El 


a. Solution 

using 

Kg 

and 

Kg 






12 

6L 

-12 

6L 


El 


6L 

4L 2 

-6L 

CM 

-3 

CM 

[Ke] 

= 7 







L 3 


-12 

-6L 

12 

-6L 


2 2 

2L^ -6L Au 


6L 
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Figure 15 Simply-Supported Beam Buckling Problem 
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[Kg]* 


30L 


36 

3L 

-36 

3L 

3L 

4L 2 

-3L 

-L 2 

-36 

-3L 

36 

-3L 

3L 

-L 2 

-3L 

4L 2 


The active degrees of freedom provide the following matrix 
equation 


El 

4 2 

Z 1 

2 4 


- P 


30L 


4 

•1 


-1 

4 




0lL 


0 



0 2 L 


0 


PL' 


Let P = 


30EI 


Then 


(4 - 4P) (2 + P) 

/% A 

(2 + P) (4 - 4 P) 


= 0 


(4 - 4P ) ( 4 - 4P) - ( 2 + PT = 0 


16 - 32P + 16P 2 - 4 - 4P - P 2 = 0 


15P* - 36P + 12 = 0 


5P 2 - 


12+ 4 = 0 



Ill 


P = 0.4, 2 

Therefore, 

(0 . 4) (30)EI 



Error = ^ = 21.6 % 

Consider a two-element model of the beam in Figure 15. 

4L 2 -6L 2L 2 0 

-6L 24 0 6L 

2L 2 0 8L 2 2L 2 

0 6L 2L 2 4L 2 

2 2 
4L -3L -L 0 

-3L 72 0 3L 

-L 2 0 8L 2 -L 2 

0 3L -L 2 4L 2 




Active DOF: 


2, 3, 4, 6 
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The active degrees of freedom provide the following matrix 
equation : 



Letting P = Yields 

30EI 



8,100 P 4 - 28,800 P 3 + 29,664 P 2 - 9,216 P + 576 = 0 
P = 0.0828653, 0.4, 1.07269, 2 


P 


30 El 
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30EI - 
P = P 

L 2 

A 

But L = L/2 

Thus 

(30) El (0.0828653) (4) 



P = 9.9438 

x 2 - 9.9438 

Error = ^ 

% 

= 0.75% 


Now consider 1 element with Ke and Kb ow . 



( 4-P) 2 - ( 2+P) 2 = 0 
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16-8P+P 2 -4-4P-P 2 = 0 

A 

12-12P = 0 
* 

P = 1 


Therefore, 

P = 12 El Which is identical with the one element 

2 

L . solution using Ke and Kg. 


Now consider the two-element model. 
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Active DOF: 2, 3, 4, 6 

The active degrees of freedom provide the following matrix 
equation : 



Letting P = yields 

12 El 

576 P 2 - 1,152 P + 576 = 0 

A 

P = 1 



12 El 
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P 


12 El P 


L 


2 


But L = L/2 


Thus 


P = 


(12) El (1) ( 


4) 


48 El 

, no good. 

It should be noted that the problem actually modeled 
by this 2 element formulation is as shown in Figure 16, 
which does not correctly model a directed force between the 
support nodes 1 and 3. 





CHAPTER 10 


GLOBAL FORMULATION OF BOW-STRING 

Examination of [Kb ow ] indicates that rigid body 
rotation capabilities occur due to the corresponding rows in 
the matrix relating to the end shears having all zero 
coefficients. Hence, any combination of [K^ow] and [Kg] 
will not have rigid body rotation capability. It has also 
been shown that assemblage of [K^ow^ elements did not model 
a force directed between the end nodes. 

Examination of a 2-element model (Appendix A) showed 
that the only fictitious forces that occurred during rigid 
body rotation were the end shears required for equilibrium. 
The shear at the center node was zero. The corresponding 
row in the geometric stiffness matrix is full, indicating 
that there is a relationship between the stiffness terms at 
each degree of freedom, and the shear at the center node. 

Examination of the first and fifth rows, however, 
would indicate that there is not any relationship between 
the end nodes. This is inherent in the assembly process, 
and is contrary to the basic supposition that we are 
considering a problem where the applied forces remain 
directed between the end nodes. 
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2 


1 



4 


3 



6 


5 




o 


P 


V 3 


o 



Figure 17 Directed Force 2-element Representation 
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Consider again Argyris's methodology for the 
directed force problem (Figure 17). 


V 3 - Vi 

Let 0 * 

L 


If one neglects the change in the axial component of 
P that occurs during rotation, as is customarily done (P cos 
0 * P), we obtain the consistent geometric stiffness matrix. 

Suppose we retain the vertical component, P sin 0, 
and use Argyris's approach to develop a load correction 
matrix. The load vector for this force becomes 


.DFC 


= P sin (V 3 -V 1 )/L, 0, 0, 0, -P sin (V 3 -V 1 )/L, 0 


The load correction matrix is generated using the equation 


[K 


DFC j 


dR 


DF 


du.; 


which yields 
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DFC 

For small rotation, cos (Vi~V 3 )/L % 1, and [K ] becomes 



DFC 

At this point, combine [Kg] + [K ] and check rigid 

DFC 

body rotation. Since [K ] contains non-zero terms in rows 
1 and 5, and [Kg] pseudo-forces occur only in the same two 
rows, only these two rows must be checked for rigid body 
rotation capability. 
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Row 1 

P j^l 2 / 5 L - 1/L , 1/10, -12/5L, 1/10, 1/L, oj -L8/2 

6 

0 

8 

L8/2 

8 


= P -1.28 + 0.58 + 0.18 + 0.18 + 0.58 
= 0 . 

Row 5 





CHAPTER 11 


PERFORMANCE OF [K T ] + [K] DFC 

DFC 

Consider a two-element model using [K>p] + [K] , 

utilizing the computer program NLBO.FOR. Table 13 compares 
the results from NLBO.FOR with the finite element solution 
using [K ] only (consistent [Ke] and [Kg] matrices). Note 
that the stiffness matrix generated by NLBO.FOR does possess 
the additional zero eigenvalue, required for a complete set 
of rigid body modes. The other frequencies have extremely 
close correlation with the traditional finite element 
solution obtained using NLFINITE . FOR . Most of them were 
identical. The largest difference was 2.8 % for frequency 
# 6 . 


When one considers the stiffness matrix [K>j>] 
generated by NLFINITE. FOR for this problem, the following is 
obtained : 
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Table 13 Frequency Comparison using NLFINITE . FOR 
and NLBO . FOR 


2 elements 


4 elements 


Freq NLFINITE. FOR 


2 E 1 em . 
0 
0 

1390 

3532 

7002 

7847 

14003 

17683 

27189 


48 IN 

30 x 10 6 PSI 
1000 IN 4 

0.03525 LB-SEC 2 /IN 2 
100 IN 

60,000,000 LBS 


NLBO. FOR 


4 E 1 em . 

2 E 1 em . 

0 

0 

0 

0 

1385 

0 

3524 

3532 

6514 

7002 

7163 

7657 

12565 

14003 

14003 

17683 

21782 

27090 

22755 


28006 


33395 


51048 


84645 


93157 



% Diff %Diff 


0 0 

0 0 

0 0 

3524 0 

6524 0 

6969 2.4 | 2.7 

12565 

14003 0 

21715 0.4 

22755 
28006 
33395 
51083 
84645 
93200 
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[K>j>] x [Rigid Body Modes] =00 0 

0 0 -6 x 10 7 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 6 x 10 7 

0 0 0 


As expected, large pseudo-forces occurred during rigid body 
rotation . 

The eigenvalues and eigenvectors generated by 
NLFINITE were: 

Lambda (1) = 0.0000 

Omega (1) = 0.0000 RAD/SEC 

The associated eigenvector is: 

0.1000000000D+01 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

0.1000000000D+01 

O.OOOOOOOOOOD+OO 


O.OOOOOOOOOOD+OO 



0.1000000000D+01 


O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

Lambda (2) = 0.0000 
Omega (2) = 0.0001 RAD/S 

The associated eigenvector is 

O.OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.1022543458D-17 
O.OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.110225609 3D-15 
O.OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.5808321841D-16 

Lambda (3) = 1930958.5265 
Omega (3) = 1389.5893 RAD/S 

The associated eigenvector is 

O.OOOOOOOOOOD+OO 

-0.1000000000D+01 

0.1508421072D-01 

O.OOOOOOOOOOD+OO 


0.1746464363D-14 



0.2418478835D-01 


O.OOOOOOOOOOD+OO 

0.1000000000D+01 

0.1508421072D-01 

Lambda (4) = 12477499.9590 
Omega (4) = 3532.3505 RAD/S 

The associated eigenvector is 

0.0000000000D+00 
0.1000000000D+01 
-0.4125765357D-01 
0.0000000000D+00 
-0.6561862202D+00 
-0.6954552347D-17 
O.OOOOOOOOOOD+OO 
0 . 1000000000D+01 
0. 41257 65357D-01 

Lambda (5) = 49021276.5957 
Omega (5) = 7001.5196 RAD/S 

The associated eigenvector is 

0.1000000000D+01 

O.OOOOOOOOOOD+OO 


O.OOOOOOOOOOD+OO 



0.2041110487D-15 


O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

-0.1000000000D+01 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

Lambda (6) = 61570298.8787 
Omega (6) = 7847.6744 RAD/S 

The associated eigenvector is 

O.OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.8622993042D-01 
O.OOOOOOOOOOD+OO 
0.1067276729D-15 
0.7402023912D-01 
0 . 0000000000 D+ 00 
-0.1000000000D+01 
-0.8622993042D-01 

Lambda (7) = 196085106.3830 
Omega (7) = 14003.0392 RAD/S 

The associated eigenvector is 

0.1000000000D+01 


O.OOOOOOOOOOD+OO 



O.OOOOOOOOOOD+OO 


-O.IOOOOOOOOOD+Ol 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

O.IOOOOOOOOOD+Ol 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

Lambda (8) = 312696968.0165 
Omega (8) = 17683.2397 

The associated eigenvector is 

O.OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.1694400209D+00 
O.OOOOOOOOOOD+OO 
0. 41200017 44D+00 
-0.361944324 3D -18 
O.OOOOOOOOOOD+OO 
O.IOOOOOOOOOD+Ol 
0.1694400209D+00 

Lambda (9) = 739222146.5762 
Omega (9) = 27188.6400 


The associated eigenvector is 
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0.0000000000D+00 
0.1000000000D+01 
-0.2124292332D+00 
0.0000000000D+00 
0.3363477755D-15 
-0. 1091964722 D+ 00 
O.OOOOOOOOOOD+OO 
-0.1000000000D+01 
-0.2124292332D+00 


Let [ 0 ] 
(eigenvectors) . 


be the matrix of mode shapes 
T 


Then [0] [K] [0] would yield a 
diagonalized stiffness matrix if all of the rigid body modes 
were present. Performing that matrix multiplication yields 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1.92 

-0 . 66 

0 

1.98 

0 

0 . 41 

-1 .36 

0 

1 

2.54xl0 6 

0 

0 

109 

0 

-1 

41 

0 

-2 

1 1 

. 19xl0 7 

0 

-3 

0 

220 

-2 

0 

0 

0 

0 5 

. 7 6xl0 7 

0 

-2 

0 

0 

0 

-1 

113 

-1 

0 4 

. 91xl0 7 

0 

-1 

-909 

0 

0 

0 

0 

-1 

0 2 

. 30xl0 8 

0 

0 

0 

2 

1 

219 

0 

-1 

0 

1.33x10 

8 -1 

0 

1 

43 

-1 

0 

-908 

0 

0 3. 

30xl0 8 
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It should be noted that small errors occur during 
the computations (initial data errors, roundoff errors, 
truncation errors, relative errors, etc.). The 2.54x10® 
term in the 3,3 position is due to the lack of rigid body 
rotation capability. The other non-diagonal terms should 
also be zero, but may be attributed to the above mentioned 
errors. The largest of these, ± 909, is still relatively 
insignificant compared to the magnitude of the diagonal 
terms . 


If one neglects the relatively small terms due to 
arithmetic errors, the following diagonal matrix is 
obtained : 


0 


0 


0 


0 0 


0 


0 0 


0 


0 0 
0 0 
0 0 
0 0 
0 0 
0 0 
0 0 
0 0 


0 0 0 

2.54x10® 0 0 

0 1 . 19xl0 7 0 

0 0 5 . 7 6xl0 7 


0 

0 

0 

0 


0 


0 0 4 . 91xl0 7 


0 

0 

0 

0 

0 


0 0 0 


0 2 . 30xl0 8 


0 

0 

0 

0 

0 

0 


0 0 0 


0 0 1 . 33x10® 


0 

0 

0 

0 

0 

0 

0 


0 0 0 


0 0 


0 3 . 30x10® 


Now consider the modified finite element solution 
from NLBO.FOR, which utilized the directed force correction 
matrix. The two-element stiffness matrix generated is: 
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. 288E8 0 0 - . 288E8 00000 

0 . 372E7 . 78E8 0 -.432E7 .78E8 000 

0 . 78E8 . 28E10 0 -.78E8 .11E10 000 

- . 288E8 0 0 .576E8 0 0 -.288E8 0 0 

0 - . 432E7 - . 78E8 0 .864E7 0 0 -.432E7 .78E8 

0 . 7 8E8 . 11E10 0 0 .56E10 0 -.78E8 .11E10 

0 00 - . 288E8 0 0 .288E8 0 0 

0 . 6E6 0 0 - . 432E7 -.78E8 0 .372E7 -.78E8 

0 000 . 78E8 . 11E10 0 -.78E8 .28E1 


The rigid body rotation matrix is 


1 


0 


0 10 0 


1 0 



T 


0 1 0010010 
0 -50 1 0 0 1 0 50 1 


[ K<p ] x [Rigid Body Modes] 


0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
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7 

The large, erroneous term (± 6 x 10 ) due to lack of rigid 
body rotation capability has been eliminated. 

The eigenvalues and eigenvectors generated by NLBO, 
corresponding to all rigid body and elastic modes and 
frequencies, were: 

Lambda (1) = -0.0122 
Omega (1) = 0.0000 RAD/S 

The associated eigenvector is: 

O.OOOOOOOOOOD+OO 
-0.9999999660D-01 
0.0000000000D+00 
0. 17017 62335D-07 
0.1999999966D-01 
0 . 0000000000D+00 
0.1000000000D+01 
0.1999999971D-01 

Lambda (2) = 0.0000 
Omega (2) = 0.0000 RAD/S 

The associated eigenvector is: 

0.0000000000D+00 

0.1000000000D+01 

-0.1021057361D-08 


0.0000000000D+00 



134 


0.9999999489D+00 
-0. 1021057 341D-08 
O.OOOOOOOOOOD+OO 
0.999999897 9D+00 
-0.1021057345D-08 

Lambda (3) = 0.0000 
Omega (3) = 0.0001 RAD/S 

The associated eigenvector is: 

0.1000000000D+01 

O.OOOOOOOOOOD+OO 

0.0000000000D+00 

0.1000000000D+01 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

0.1000000000D+01 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

Lambda (4) = 12477499.9590 
Omega (4) = 3532.3505 RAD/S 

The associated eigenvector is: 

O.OOOOOOOOOOD+OO 


0.1000000000D+01 
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-0. 41257 65357D-01 
0.0000000000D+00 
-0.6561862202D+00 
-0.13734 47 533D-16 
0.0000000000D+00 
0.1000000000D+01 
0. 41257 65357D-01 

Lambda (5) = 49021276.5957 
Omega (5) = 7001.5196 RAD/S 

The associated eigenvector is: 


0.1000000000D+01 

0.0000000000D+00 

O.OOOOOOOOOOD+OO 

0.8376574057D-16 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

-0.1000000000D+01 

O.OOOOOOOOOOD+OO 

O.OOOOOOOOOOD+OO 

Lambda (6) = 58630070.5158 

Omega (6) = 7657.0275 

The associated eigenvector is: 


O.OOOOOOOOOOD+OO 



0.1000000000D+01 


-0. 895141145 4D-01 
0.0000000000D+00 
0.125912 67 9 4D -15 
0.7572882821D-01 
0.0000000000D+00 
-0.1000000000D+01 
-0. 895141145 4D -01 

Lambda (7) = 196085106.3830 
Omega (7) = 14003.0392 RAD/S 

The associated eigenvector is 

0.1000000000D+01 
0.0000000000D+00 
0.0000000000D+00 
-0.1000000000D+01 
0 .0000000000D+00 
0 .0000000000D+00 
0.1000000000D+01 
0 .0000000000D+00 
0.0000000000D+00 

Lambda (8) = 312696968.0165 
Omega (8) = 17683.2397 RAD/S 


The associated eigenvector is 
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O.OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.1694400209D+00 
O.OOOOOOOOOOD+OO 
0.41200017 4 4D+00 
0. 328045148 6D -16 
O.OOOOOOOOOOD+OO 
0.1000000000D+01 
0.1694400209D+00 


Lambda (9) = 733880567.5204 
Omega (9) = 27090.2301 RAD/S 

The associated eigenvector is: 

O.OOOOOOOOOOD+OO 

0.1000000000D+01 

-0.2134858855D+00 

O.OOOOOOOOOOD+OO 

0.6206811895D-15 

-0.1102288283D+00 

O.OOOOOOOOOOD+OO 

-0.1000000000D+01 

-0.2134858855D+00 


Let [ 0 ] 


be the matrix of mode 


( eigenvectors ) . 


[0] T [K] [0] should 


shapes 
i e 1 d a 


Then 
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diagonalized stiffness matrix if all of the rigid body modes 


were 

present . 

Performing 

that 

matrix 

multiplication 

yi el ds 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0.08 

0 

0 

0 

0.34 

0 

0 

1.71 

1 

-0.55 

0 

-1.97 

0 

0.33 

0.84 

0 

1 

-1.99 

1.19E7 

0 

-3 

0 

-39 

-1 

0 

0 

0 

0 5. 

7 6E7 

0 

-2 

0 

0 

0 

-1 

-1.91 

-2 

0 

4. 91E7 

0 

-1 

-81 

0 

0 

0 

0 

-1 

0 

2.31E8 

0 

0 

0 

0 

0 . 56 

-39 

0 

1 

0 

1.33E8 

-1 

0 

0 

1.65 

-1 

0 

-81 

0 

-2 

3. 32E8 


It should be noted that minor errors still occur 
during the computations (initial data errors, roundoff 
errors, truncation errors, relative errors, etc.)- The 
examination of these errors is beyond the scope of this 
dissertation. It can be readily seen, however, that the 
largest of these error has been reduced an order of 
magnitude (from ± 909 to 81). 

Neglecting the relatively small terms due to 
arithmetic errors, the following diagonal matrix is 


obtained : 
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0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

1.19E7 

0 

0 

0 

0 

0 


0 

0 

0 

0 

5.76E7 

0 

0 

0 

0 


0 

0 

0 

0 

0 

4.91E7 

0 

0 

0 


0 

0 

0 

0 

0 

0 

2.31E8 

0 

0 


0 

0 

0 

0 

0 

0 

0 

1.33E8 


0 

0 

0 

0 

0 

0 

0 

0 


3.32E8 


Most importantly, the large erroneous term in the 
3,3 position of the matrix obtained using the conventional 
finite element formulation is now identically zero, and the 
matrix has been properly diagonalized. Thus, adding [K] , 
as developed in this dissertation, corrected the lack of 
rigid body rotation capability of the pre-loaded beam 
element, as well as provided the correct diagonalized 
stiffness matrix in the d i agona 1 i z a t i on/ pa r t i t i oni ng 
methodology used in finite element dynamic analysis. 



CHAPTER 12 


SUMMARY 

Based upon this investigation, the following 

conclusions have been developed: 

1. Grounding is due to the development of pseudo-forces at 
the element level required to counteract a force- 
imbalance inherent in the development. This causes a 
lack of rigid body rotational capability of the 
geometric stiffness matrix. 

2. Although the consistent geometric stiffness matrix 
provides acceptable results for most static displacement 
and buckling problems, provided a sufficient mesh is 
used, modifications of the global stiffness matrix 
(zeroing out of erroneous terms, and appending the 
missing rigid body modes) must be done to more 
accurately predict the dynamic response. 

3. Although the rigid body mode test is routinely used to 
detect the presence of modeling errors in finite element 
models, it is not sufficient reason to invalidate a 
model subjected to pre-loads. 
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4. Various higher order stiffness matrices developed by 
others, which include shear and rotatory inertia 
effects, were examined. As expected, the inclusion of 
these higher order effects does not compensate for the 
inaccuracy (lack of rigid body rotation capability) of 
the geometric stiffness matrix. 

5. Rigorous solutions of the pre-loaded beam with various 
end conditions were developed. 

6. The Galerkin criterion was used to develop stiffness and 
mass matrices from the rigorous solutions, which were 
incorporated into a modified finite element algorithm. 

7. Sample problems involving pre-loaded beams with various 
spring support conditions were solved using the modified 
finite element algorithms, and the results compared with 
the rigorous solutions. 

8. The occurrence of dynamic "flutter" instabilities was 
determined by the rigorous solutions. There was good 
correlation obtained using the modified finite element 
algorithms . 

9. The tangential stiffness matrices developed did not 
possess the three zero eigenvalues required for all the 
rigid body modes. 
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10. The directed force (ie., bow-string) problem was 
examined, since the force unbalance inherent in the 
other developments does not occur in this situation. 

11. Development of the bow-string stiffness using Clough and 
Penzien's technique provided a modified matrix, but it 
was shown to lack rigid body rotational capability. 

12. Development of the bow-string stiffness using Saunders 
methodology provided a modified matrix, but it also was 
shown to lack rigid body rotational capability. 

13. A bow-string stiffness matrix developed from the 
rigorous solution using Galerkin's criterion possessed 
all the required rigid body modes. [Kgow3 was shown to 
provide an acceptable first approximation to the 
directed force buckling problem solved by Timoshenko and 
Gere, and it performs properly in the 
diagonal ization/partitioning methodology used in dynamic 
response. However, it does not properly model a pre- 
loaded beam where the force is directed between the end 
nodes only of an assembled mesh. 

14. By considering the directed force problem at the global 
level, using traditional development of [Kg] from the 
horizontal component of the directed force, and 
Argyris's load correction method for the vertical 
component, a load correction matrix [K ] was 
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developed, which, when combined with [Kg], provided a 
complete set of rigid body modes. This combined matrix 
performs properly in the diagonal ization/partitioning 
methodology used in dynamic response. 

There is the potential for a great deal of future 

work with the directed force beam element and the technique 

DFG 

used in its development. The use of [Kg] + [K ] should be 
compared with the results using Craig-Bampton ' s 
substructuring scheme for various beams. In addition, 
physical testing of a pre-loaded directed force beam with 
free/free boundary conditions should be undertaken for 
comparison. The incorporation of these techniques in the 
development of a directed force correction matrix for pre- 
loaded membrane elements would be a logical extension. 
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ABSTRACT 

Deployable solar arrays consist of a "blanket" of 
solar collectors, and a mast. The blanket is stretched into 
position when the array is deployed. The stiffness of the 
array is a function of the rigidity of the mast as well as 
the tension in the blanket. 

Current finite element frequency analysis consists 
of using MSC Nastran solution 64 (non-linear analysis) to 
obtain the tangential stiffness matrix of the array. This 
matrix is then input, using DMAP alters, into MSC/Nastran 
Solution 63 (dynamic analysis) to obtain the natural 
frequencies of the array. 

The author has found that ps eudo - f o r c es are 
developed, however, at the element level due to limitations 
inherent in the geometric stiffness matrices currently in 
accepted use. In particular, the geometric stiffness 
matrices lack the capability for rigid body rotation, 
especially when the rotations are large. 

The author demonstrates the limitations of the 
analysis, shows where the errors are introduced in the 
derivation of the geometric stiffness matrix, and examines 
various techniques either to eliminate the pseudo-force 
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generation and/or improve upon the convergence of the 
current algorithms. 

* 

This paper is the product of an NASA/ASEE Summer 
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between Cleveland State University and the NASA Lewis 
Research Center. 


* 


National Aeronautics and Space Administration/American 
Society of Engineering Educators 



152 


SPACE STATION SOLAR ARRAY 

NASA's space station is powered utilizing 
photovoltaic arrays, consisting of a deployable truss "mast" 
and blanket substrates (Figure 18). The stiffness of the 
split-blanket array is a function of the rigidity of the 
mast as well as the tension maintained in the blankets. The 
"blankets" themself possess negligible stiffness. 

The free vibration characteristics of the split- 
blanket solar arrays was studied using two methods. Mode 
shapes and frequencies were calculated using equations of 
continuum mechanics, as well as a finite element solution 
using MSC Nastran [1] and [2]. 

The finite element modeling consisted of generating 
a tangential stiffness matrix by applying the pre-tensioning 
load in MSC/Nastran geometric non-linear solution (solution 
64). The stiffness matrix generated was then input into 
MSC/Nastran dynamic analysis (solution 63) to obtain the 
natural frequencies and mode shapes [3]. 

The finite element analysis indicated that large 
internal "pseudo-forces" developed when rigid body motion 
was applied. An investigation was subsequently made to 
determine whether the large pseudo-forces which developed 
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Figure 18 Space Station Split-Blanket Solar Array 
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were due to user modeling errors or limitations of the 
finite element process. 

The geometric stiffness matrix utilized in 
MSC/Nastran solution 64 was found to be identical to the 
stiffness matrix formulated by Martin [4], For simplicity, 
a 2-dimensional beam-column element was investigated. 

LIMITATIONS OF CURRENT [Kg] 

Typically, finite element static analysis is used to 
solve linear elastic problems of the form 

[Ke] {u} = {R} --- - (1) 

where [Ke] is the elastic stiffness matrix 
{u} is the nodal displacement vector 
{ R } is the force vector 

The [Ke] matrix must possess the capacity for rigid 
body displacement. In other words, the element must be able 
to translate or rotate without developing stresses (Figure 
19). 


The rigid body translation vector is {v, 0, v, 0). 
Similarly, the rigid body rotation vector is approximated by 
{-L0/2,0,L0/2,0} , which can be written as 0{ -L/2 , 1 , L/2 , 1} , 
where 0 is the angle of rotation. Combination yields 
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1 3 



a. 2-Node Bernoulli Beam Element with 
Degrees of Freedom Shown 




Figure 19 Rigid Body Modes 
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rigid 

body 

modes 


v -L9/2 
0 8 

v L9/2 
0 0 


The elastic stiffness matrix for a 2-node Bernoulli 

beam is 


[Ke] 


12 

6L 

-12 

6L 

6L 

4L 2 

-6L 

CM 

iJ 

CM 

-12 

-6L 

12 

-6L 

6L 

2L 2 

-6L 

CM 


where E = modulus of elasticity 
I = moment of inertia 
L = length of element 


( 2 ) 


By definition of rigid body motion, 

[Ke} {rigid body modes) must equal {0}. (3) 

Substituting (1) and (2) into (3) yields 


12 

6L 

-12 

6L 

6L 

4L 2 

-6L 

2L 2 

-12 

-6L 

12 

-6L 

6L 

2L 2 

-6L 

CM 

1-3 


1 -L/2 


0 0 

0 1 


0 0 

1 L/2 


0 0 

o 

L. 


— — - , 
o 

O 

i 
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Therefore, [Ke] has the capacity for rigid body motion. 

NASA LeRC routinely uses a DMAP alter which 
calculates 

[K] {rigid body modes} = {RFORCES (pseudo-forces)}. 

For an elastic problem, the presence of large RFORCES would 
indicate that stresses are being produced during rigid body 
movement. These pseudo-forces are an indication that 
"grounding" has occurred, and that the model is not 
reliable. 

Many problems, such as the solar array, are non- 
linear problems. Finite element solves non-linear problems 
of the form 

[[Ke] + [Kg]] {u} = {R} - {F} 

where [Ke] is the elastic stiffness matrix. 

[Kg] is the geometric stiffness matrix. 

{R} is the output force vector at the end of a step. 
{F} is the input force vector at the beginning of a 
step . 

{u} is the change in the displacement vector during a 
step . 

Graphically, this is shown by Figure 20. 

The traditional [Kg] matrix, developed by Martin 


[4], is 
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6/5L 1/10 -6/5L 1/10 

1/10 2L/15 -1/10 -L/30 

[Kg] = P 0 

-6/5L -1/10 6/5L -1/10 

1/10 -L/30 -1/10 2L/15 

If I apply rigid body displacement during an incremental 
load step, I obtain 

6/5L 1/10 -6/5L 

1/10 2L/15 -1/10 

Po 

-6/5L -1/10 6/5L 

1/10 -L/30 -1/10 

It can be seen that [Kg] possesses the capacity for 
rigid body translation, but not rigid body rotation. Thus, 
Martin's [Kg] is not exact. It can also be shown that 
MSC/Nastran non-linear analysis (based on Martin's 
development) similarly does not have an exact geometric 
stiffness matrix and will produce pseudo-forces. Therefore, 
the RFORCE check DMAP alter for MSC/Nastran solution 64 
(non-linear) analysis) is not sufficient criteria for 
determining the validity of a model. 

In spite of its deficiencies, Martin's [Kg] provides 
acceptable results for solving statics problems due to the 
iteration process. (Although [Kg] is not exact, the process 
converges to the exact solution.) 


1/10 1 -L8/2 0 -P o 0 

-L/30 0 0 0 0 

-1/10 1 L0/2 0 P o 0 

2L/15 0 0 0 0 
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In dynamic analysis, however, finite element is 
used to solve equations of the form 

[M] {u’} + [K] {u} = {R} 

where [M] is the mass matrix. 

{u} is the displacement (mode shape) vector. 

{u } is the second derivative with respect to time of 
the displacement vector. 

[K] is the stiffness matrix (Ke or Ke + Kg, when 
applicable) . 

(R) is the excitation forces 

There is no apparent guarantee that the natural frequencies 
of vibration from (4) are accurate, when [Kg] is known to be 
inexact . 


LARGE ROTATION EFFECTS 

The rigid body rotation vector previously used is 
C { - L 9 / 2 , 0 , L 0 / 2 , 9 } . Figure 21 illustrates rigid body 
rotation of a beam with an axial load [5]. If we let the 

A ’ 

angle of rotation equal 2fl, then u = B{-L,2,L,2} T . 

Consider the work/energy relationship of Figure 21 
work done by P G = P 0 L( 1-cos 2fl) 

= 2P 0 L(l-cos 2fl)/2 

But, (1-cos 2fl)/2 = sin 2 B * B 2 + 0(B 4 ) + ••• 

A A 

Therefore, the work done = 2P 0 B 2 + 0(fl 4 ) = -V = u T K u/2, 
where -V is the loss of potential energy. 




Figure 21 Rigid Body Rotation Angle of 2fl 
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Similarly 


u T K u/2 = P 0 J3 2 / 2 [ -2 0 2 0] 


■L 

2 


P 0 B 2 (4L) = 2P 0 LB 2 


L 

2 


Therefore, Martin's [Kg] provides a correct energy 
relationship for the representation shown. 


It should be noted, however, that the displacement 
in the Y direction (axial direction in the original 
geometry) has been neglected when we let 



which is really 


0 

-L 

2 

0 

L 

2 


(5) 


The zero terms in (5) negate any contribution to the 
equation from axial terms in the stiffness matrix. (If 
there were any axial terms.) Since significant axial 
loading occurs in the solar array (as well as other 
beam/column problems), and these axial loads significantly 
affect the stiffness, it intuitively seems unreasonable to 
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loading occurs in the solar array (as well as other 
beam/column problems), and these axial loads significantly 
affect the stiffness, it intuitively seems unreasonable to 
arbitrarily neglect the contribution of axial stiffness 
terms and axial displacements. 

The exact rigid body rotation vector is 


u exact 


L/2 (1-COS 213) 
-L/2 SIN 2B 
20 

-L/2 (1-COS 20) 
L/2 SIN 20 
20 


( 6 ) 


Series expansion, and truncation of higher order terms 
yields 


u = 


L0“ 

-Lfl 

20 

-Lfl" 

L0 


(7) 


20 
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0 0 0 0 

0 6/5L 1/10 0 

0 1/10 2L/15 0 

Po 

0 0 0 0 

0 -6/5L -1/10 0 

0 1/10 -L/30 0 

and 

AL 2 /I 0 0 -AL 2 /I 0 0 

0 12 6L 0 -12 6L 

0 6L 4L 2 0 -6L 2L 2 

-AL 2 /I 0 0 AL 2 /I 0 0 

0 -12 -6L 0 12 -6L 

0 6L 2L * 0 -6L 4L* 

From (8) and (9), when the more exact rigid body 
rotation vector is used, neither [Ke] nor [Kg] possesses 
rigid body rotation capabilities, although from (8), 2AEQ 2 
approaches 0 as half the angle of rotation gets very small. 
A similar procedure shows that both possess rigid body 
translation capability in two directions. 
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DEVELOPMENT OF MODIFIED [Kg] 

Using rigid body translation relationships, and the 
expanded rigid body rotation vector, I can solve for 
additional terms in the [Kg] matrix which enforces the rigid 
body capabilities. 

Using rigid body translation constraints, [Kg] has 

the form 


[Kg_l=P c 


A 

B 

C 

-A 

-B 

C 

B 

6/5L 

1/10 

-B 

-6/5L 

1/10 

C 

1/10 

2L/15 

-c 

-1/10 

-L/30 

-A 

-B 

-C 

A 

B 

-C 

-B 

-6/5L 

-1/10 

B 

6/5L 

-1/10 

C 

1/10 

-L/30 

-C 

-1/10 

2L/15 


( 10 ) 


Multiplying (10) by (7) and setting the product 
equal to the zero vector yields 


A 

B 

C 

B 

6/5L 

1/10 

C 

1/10 

2L/15 

-A 

-B 

-C 

-B 

-6/5L 

-1/10 

C 

1/10 

-L/30 


-A -B C 

-B -6/5L 1/10 

-C -1/10 -L/30 

A B -C 

B 6/5L -1/10 

-C -1/10 2L/15 



LB 2 


0 


-Po 


0 


2B 


0 


-LB 2 


0 


BL 


0 


2B 


0 


(ID 


Expanding row 3 yields 
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P 0 [2CLB 2 - BL/5 + 4BL/15 - 2BL/30] = 0 

P o 2CL0 2 = 0 
C = 0 


Expanding row 2 yields 

P Q [ 2BLG 2 - 12BL/5 + 4A/10] = 0 

2BLB -2=0 
B = 1/Lfl 


Expanding row 1 yields 

P Q [2ALA 2 - 2] = 0 
ALA 2 -1=0 
A = 1/LB 2 


Thus, substituting (12), (13), and (14) into (10) yields 


[KslL=P 0 


1/Lfl 2 

1/Lfl 

0 

-1/LB 2 

-1/Lfl 

0 

1/LB 

6/5L 

1/10 

-1/Lfl 

-6/5L 

1/10 

0 

1/10 

2L/15 

0 

-1/10 

-L/30 

-1/LB 2 

-1/Lfl 

0 

1/Lfl 2 

1/Lfl 

0 

-1/Lfl 

-6/5L 

-1/10 

1/Lfl 

6/5L 

-1/10 

0 

1/10 

-L/30 

0 

-1/10 

2L/15 


( 12 ) 


(13) 


(14) 


(15) 



167 


The strain energy U should equal zero. It can be 
calculated from the equation 

P Q fl 2 / 2 [LB,-L,2,-LB,L,2] [Kg.] LB 

-L 

2 

(16) 

-LB 

L 

0 

2 

Performing the matrix multiplication in (16) yields 
T 

[0, 0, 0, 0, 0, 0] . Therefore, no strain energy occurs 

during rigid body rotation using [Kg.] . 
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AL 2 /I k 12 K 13 -AL 2 / I -K 12 K 16 0 

K 12 12 6L -K 12 -12 6L 0 

El K 13 6L 4L 2 -K 13 -6L 2L 2 0 

3 2 2 = (1®) 

L° -AL / I -K 12 -K 13 AL^/I K 12 -K 16 0 

-K 12 -12 -6L K 12 12 -6L 0 

K 16 6L 2L* -K 16 -6L 4L* 0 

Expanding row 2 yields 

2K 12 L0 2 - 24LB + 24LB = 0 

2K 12 LB 2 = 0 (19) 

Therefore, K^ 2 must equal 0. 

Expanding row 3 yields 

-2K 13 Lfl 2 - 12L 2 B + 8L 2 B + 4L 2 fl = 0 

-2K 13 L0 2 = 0 (20) 

Therefore, K^ 3 must equal zero. 

If Ki 2 and K ]_ 3 equal zero, no additional 

coefficients appear in line 1 of the [Ke] matrix. Suppose I 

add a correction term to K^. 

(AL 2 /I+C)Lfl 2 -(AL 2 /I+C)(-Lfl 2 ) = 2 ( AL 2 / I+C) ( LB 2 ) (21) 

From (21) it can be seen that adding a correction term to 
Kn wi 11 not eliminate the error. 

Thus, [Ke] can not be modified to obtain rigid body 
rotation capabilities for large rotations utilizing the 
procedure used to modify [Kg]. The only possibility for 
improving [Ke] must include corrections to existing terms. 
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VERIFICATION OF MODIFIED [Kg] 

Several tests of the modified stiffness matrix were 
undertaken. Summary of the tests follow, 
a. RFORCE Check 

One and two element beam stiffness matrices were 
multiplied with the expanded rigid body matrix. In all 
cases, no pseudo-forces were produced from rigid body 
translation. The expanded rigid body rotation vector 
yielded the following pseudo-forces. 

Matrix Pseudo-forces 

[Ke] one element [2Afl 2 E,0,0,-2Afi 2 E,0,0] 

two elements [2Afl 2 E,0,0,0,0,0,0,-2Afl 2 E,0,0] 

[Kg] one element [ 0 , -2BP 0 , 0 ,0 , 2BP 0 , 0 ] 

two elements [0 , -2BP 0 , 0 ,0 ,0 ,0 ,0 , 2BP 0 ,0 ] 

[Kg.] one element [0,0,0,0,0,03 

[Kg.] two elements [0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] 

[Ke+Kg] one element [ 2AB 2 E , -2 BPq , 0 , -2AB 2 E , 2BP 0 , 0 ] 

[Ke+Kg] two elements [ 2AB 2 E , -2BP 0 , 0 ,0 , 0 , 0 , -2AB 2 E , 2BP 0 , 0 ] 

[Ke+Kg] one element [ 2AB 2 E , 0 , 0 , -2Afl 2 E, 0 , 0 ] 

[Ke+Kg] two elements [ 2AB 2 E , 0 , 0 , 0 , 0 , 0 , -2AB 2 E , 0 , 0 ] 

Based on the above, the modified [Kg] matrix 

eliminates the fl error terms generated using [Kg] standard. 

2 

The B terms generated from [Ke] remain. Thus, when B (half 
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the angle of rotation) is less than one radian, the total 
error is reduced. The error is associated with a pseudo- 
axial force only. 

b. Eigenvalues 

Using [Ke] and [Kg] standard, it was found that only 
two zero eigenvalues exist. [ILSL] has three zero 
eigenvalues, plus one eigenvalues which is very close to 
zero. The eigenvalues are close to zero even when only one 
or two elements are used. 

c. Stability Analysis 

[Kg.] was used in the solution of a simply supported 
beam subjected to an axial load. The critical buckling load 
was calculated, and compared with the traditional solution 
using Martin's [Kg], as well as the exact solution. 

When the boundary condition, B^ equals -B2 ( as 
occurs during buckling) was applied, the solution using 
[Kg.] was identical with the solution using Martin's [Kg]. 

CONCLUSIONS 

Continuing effort is being made on improving the 
capabilities of the element stiffness matrices used in 
the solar array dynamic analysis. The modified [Kg] 
developed reduces the pseudo-forces produced in the beam- 
column stiffening problem, provided that the angle of 
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rotation is less than two radians. Modifications roust be 
developed, however, to eliminate the pseudo-force 
contributions from [Ke] which result from utilizing the 
expanded rigid body rotation vector. This would permit the 
tangential stiffness matrix [K-p] to possess the three zero 
eigenvalues associated with rigid body movement of any 
magnitude . 

It was disappointing that the [Kg] developed did not 
approve upon the relatively slow convergence rate of the 
stability problem. Further investigation is needed to 
determine whether a modified [Ke] + [Kg] would improve upon 
this convergence rate. 

Extension of the modifications to the stiffness 
matrices of other elements, particularly plate elements, 
will also be developed. 

Finally, testing of the performance of the modified 
matrices in the actual solar array model, and comparison 
with the continuum mechanics approach, will be performed. 
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APPENDIX B 

Matrix Methods (Dr. Bellini, Cleveland State University) 
Finite Element Approach 
% = U - V 
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Which are the shape functions for the static beam. 

Therefore, [Kg] is only approximate. 


Hi = 1 - x/L 

H 2 = (l-3(x/L ) 2 + 2(x/L) 3 ) 
H 3 = (x - 2 x 2 /L + x 3 /L 2 ) 


H 4 = x/L 

H 5 = (3(x/L ) 2 - 2 ( x/L ) 3 ) 
H 6 = (~x 2 /L + x 3 /L 2 ) 


Shape functions for the static beam (Figure 22). 


u = ui + H 4 u 4 


du 

— = Hi' ui + H 4 ’ u 4 
dx 


W = H 2 u 2 + H 3 u 3 + H 5 U5 + H 6 u 6 


dW 

— = H 2 'u 2 + H 3 'u 3 + H 5 ' U 5 + Hg'ug 

dx 


— =• = H 2 "u 2 + H 3 "u 3 + H 5 "u 5 + H 6 ”u 6 

dx 


Term (1> 
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EA 


du D 
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L 






dx 


0 


0 
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Figure 22 2-Node Element Degrees of Freedom 
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For constant EA 



H 2 " = - 6 /L 2 + 12x/L 3 
H 3 " = -4/L + 6 x/L 2 
H 5 " = 6 /L 2 - 12x/L 3 
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H 6 " = - 2 /L + 6x/L 2 
For El = constant 


1 r 
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CM 
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C-3 
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duo 

Setting EA = P 

dx 


Axial force in member, and assuming 
P = constant + tension 

- compression 


1 r 



] p 

36 

3L 

-36 

3L 
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2 L 

u 3 

u 5 

^6 — 

J 30L 
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4L 2 
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APPENDIX C 


NLFINITE.FOR Computer Program and Output 


C NLFINITE.FOR (GEOMETRIC NONLINEAR PROBS ) 

C REVISED 10-24-90 ( IF STATEMENT IN JCBI REVISED 
C MODIFICATION OF PROGRAM FINITEL .FOR VIBRATION 
C ANALYSIS OF BEAMS, RODS, AND PLANE FRAMES USING BEAM 
C ELEMENTS WITH AN AXIAL PRETENSION LOAD. SUBROUTINES 
C JCBI, DECOMP, MATINV , MATMPY , AND SEARCH. 

C DEVICE * IN READ AND WRITE STATEMENTS IS THE CONSOLE. 

C DEVICE 2 IN WRITE STATEMENTS IS THE PRINTER. 

C * IN THE PLACE OF A FORMAT STATEMENT NUMBER MEANS FREE FORMAT. 
IMPLICIT REAL*8 (A-H.O-Z) 

REAL*8 L.IA.KEL.MEL.KEG.MEG 
INTEGER SUB .ROWSUB .COLSUB ,B ,Z ,EN ,CFIX 

DIMENSION KEL( 6,6,8) ,MEL( 6,6,8 ) ,KEG( 6 ,6 .8 ) ,MEG(6 ,6 ,8 ) , 
*RT(6,6). 

*R( 6 ,6 ) ,TK( 6 ,6 ) ,TM( 6 ,6 ) ,SK( 27 ,27 ) ,SM( 27 ,27 ) ,RSK( 27,27 ) , 

*RSM( 27 ,27 ) ,E( 8 ) , A( 8 ) ,X( 9 ) , Y( 9 ) ,GAMMA( 8 ) ,IA( 8 ) 

DIMENSION JNM( 8 ,2 ) ,CFIX( 27 ) ,SUB( 6 ) 

OPEN( UNIT=2 .FILE* ’PRN * ) 

C READ IN PROBLEM DATA AS INDICATED BY MESSAGES ON CONSOLE. DATA 
C READ IN IS PRINTED OUT. (PROGRAM STATEMENTS 2 THROUGH 40) 
WRITE( 2,1) 

1 FORMAT( ’NLFINITE.FOR; FULL KE+KG MATRICES .REVISED 10-24-90') 
WRITE( * ,2 ) 

2 FORMAT( / , ’ ENTER THE NUMBER OF BEAM ELEMENTS (II)’,/) 
READ(*,3)NUMEL 

3 FORMAT( II ) 

DO 4 1*1 ,NUMEL 
• WRITE( * ,5 )I 

4 READ( *,6)A( I ) 

5 FORMAT( / , ’ ENTER A(',I1,’) (F20.0)*,/) 

6 FORMAT( F20 .0 ) 

WRITE( 2,7) 

7 F0RMAT(6X, 'THE AREA ARRAY A IS:’,/) 

DO 8 1*1 .NUMEL 

8 WRITE( 2 ,9 )I ,A( I ) 

9 F0RMAT(6X,’A( ’ ,11,’) * ’.E14.7) 

DO 10 1*1 .NUMEL 

WRITE(*,11 )I 

10 READ( * ,12 )E( I) 

11 FORMAT( / , ’ ENTER E(’,I1,’) (F20.0)',/) 

12 FORMAT( F20 .0 ) 

WRITE( 2,13) 

13 F ORMAT ( /6X , ’ THE ELASTICITY ARRAY E IS;',/) 

DO 14 1*1 .NUMEL 

14 WRITE( 2,15)1 ,E( I ) 

15 F0RMAT(6X, ’E( ' ,11 , ’ ) * ’.E14.7) 

DO 16 1 = 1 .NUMEL 
'WRITE(*,17)I 
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16 READ( * ,18 )IA( I ) 

17 FORMAT( / , 'ENTER IA( ’ ,11 , ' ) (F20.0)’ ./) 

18 FORMATC F20.0) 

UIRITE( 2 ,19 ) 

19 FORMAT( / ,6X , 'THE MOMENT OF INERTIA ARRAY IA IS;’,/) 

DO 20 1 = 1 .NUMEL 

20 URITEC 2 ,21 )I ,IA( I ) 

21 FORMATf 6X,’IA( ',11,') * ’ ,E14 .7 ) 

DO 82 1 = 1 ,NUMEL 

URITEC * ,83 )I 

82 READ( * ,84 )GAMMA( I ) 

83 FORMAT( / , ' ENTER GAMMA( ’ , 1 1 , ' ) (F20.0)'./) 

84 FORMAT( F20 .0 ) 

URITEC 2 ,85 ) 

85 FORMAT( / ,6X , 'THE GAMMA ARRAY IS;',/) 

DO 86 1=1 .NUMEL 

86 URITEC 2 ,87 )I ,GAMMA( I ) 

87 FORMAT( 6X , ’GAMMA( ' ,11 , ’ ) = ’ ,E14 .7) 

URITEC * ,88 ) 

88 FORMAT ( ’ ’, 'ENTER THE AXIAL TENSION PRELOAD (PLOAD)’,/) 

READC * ,89 )PLOAD 

89 FORMAT( F20 .0 ) 

URITEC 2 ,90 ) 

90 FORMATC / ,6X , 'THE AXIAL PRETENSION LOAD IS;’,/) 

URITEC 2, 91 )PLOAD 

91 FORMATC F20.0) 

URITEC* ,22 ) 

22 FORMATC/, 'ENTER THE NUMBER OF JOINTS, NJTS (II)’,/) 

READC *,23 )NJTS 

23 FORMATC II ) 

DO 24 1=1, NUMEL 
DO 24 J=1 ,2 
URITEC *, 25 )I,J 

24 READC *,26 )JNM( I ,J ) 

25 FORMATC/,’ ENTER JNMC ’ .11 , ’ , ’ ,11 , ’ ) (II)’,/) 

26 FORMATC II) 

URITEC 2, 27) 

27 FORMATC /,6X, ’THE JOINT-NUMBER MATRIX IS;’,/) 

DO 28, 1 = 1 .NUMEL 

28 URITEC 2 ,29 )JNM( 1,1) , JNM( 1,2) 

29 FORMATC 10X ,15 ,14 ) 

DO 30 1 = 1, NJTS 
URITEC* ,31 )1 ,1 

30 READ( * , * )X( I ) , Y( I ) 

31 FORMATC/’ ENTER JOINT COORD. X( ’ ,11 . ’ ),Y( ’ ,11 , ’ )( 2F20 .0 )’ ,/) 
URITEC 2, 33) 

33 FORMATC /,6X, ’THE JOINT COORDINATES ARE;’,/) 

DO 34 1=1, NJTS 

34 URITEC 2 ,35 )I ,X( I ) ,1 ,Y( I ) 

35 FORMATC 6X, ’X( ’ ,11, ’ )=' ,E14 .7 ,5X , ’Y( ’ ,11 , ’ )= ’ ,E14 .7 ) 

URITEC *,36) 

36 FORMATC/,’ ENTER THE NUMBER OF FIXED COORDINATES( 12 ) ’ ,/ ) 
READ( * ,37 )NB 

37 FORMATC 12) 

IF(NB.EQ.O)GO TO 94 
DO 38 1=1 ,NB 
URITEC*, 39)1 
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38 READ( * ,40 )CFIX( I ) 

39 F0RMAT( /, ’ ENTER CFIX( ’ ,12, ’ ) ( 12)’ ,/) 

40 F0RMAT( 12 ) 

WRITE( 2 ,41 ) 

41 FORMAT^ / ,6X , ’ARRAY CFIX IS;’,/) 

DO 42 1*1 ,NB 

42 WRITEC 2,43)1 ,CFIX( I ) 

43 FORMAT( 6X, ’CFIX( ’ ,12, ' )=’ ,12) 

C GENERATE NULL 3-DIMENSIONAL ARRAYS KEL AND MEL . PLANES OF KEL 
C AND MEL WILL LATER CONTAIN THE LOCAL ELEMENT STIFFNESS AND MASS 
C MATRICES, RESPECTIVELY. 

94 DO 44 1 = 1 ,6 
DO 44 J*1 ,6 
DO 44 M=1 .NUMEL 
KEL( I ,J,M)=0. 

44 MEL( I ,J,M)=0. 

C GENERATE NULL MATRICES R AND RT WHICH WILL LATER BECOME THE 
C TRANSFORMATION MATRIX AND ITS TRANSPOSE, RESPECTIVELY 
-DO 45 1 = 1 ,6 
DO 45 J=1 ,6 
R(I ,J)=0.0 

45 RT( I ,J )=0 .0 

C GENERATE THE LOCAL ELEMENT STIFFNESS MATRICES AND STORE IN THE 3 
C 3-DIMENSIONAL STIFFNESS ARRAY KEL (SEE FIG. 8-11 FOR THE 
C EQUATION USED). EACH PLANE IN THE 3-DIM. ARRAY IS ONE ELEMENT 
C STIFFNESS MATRIX. 

DO 100 EN=1 .NUMEL 
IC= JNM( EN ,1 ) 

ID=JNM( EN ,2 ) 

L=DSQRT( ( X( ID )-X( IC ) )**2+( Y( ID )-Y( IC ) )**2 ) 

QUOT=IA( EN )/A( EN ) 

R1=DSQRT( QUOT ) 

F=E(EN)*IA(EN)/L 
P=F/R1**2 
0=4 ,*P*R1**2 
S-3 .*0/( 2 . *L ) 

T=S*2./L 

SINA=( Y( ID )-Y( IC ) )/L 
COSA=( X( ID )-X( IC ) )/L 
KEL( 1,1 ,EN)=P 
KEL( 1 ,4 ,EN )=-P 

KEL( 2 ,2 ,EN )*T+( 6 .*PLOAD/( 5 .*L ) ) 

KEL( 2 ,3 ,EN )=S + . 1*PL0AD 

KEL( 2 ,5 ,EN )=-T-( 6 .*PLOAD/( 5 .*L ) ) 

KEL( 2 ,6 ,EN )»S+ . 1*PL0AD 

KEL( 3,3 ,EN )=0+( 2 ,*PL0AD*L/15 . ) 

KEL( 3 ,5 ,EN )=-S- . 1*PL0AD 

KEL( 3 ,6 ,EN )=0/2 .-( PL0AD*L/30 . ) 

KEL( 4 ,4 ,EN )=P 

KEL( 5 ,5 ,EN )=T+( 6 .*PLOAD/( 5 .*L ) ) 

KEL( 5 ,6 ,EN )=-S- . 1*PL0AD 

KEL( 6,6 , EN )=Q+( 2 .*PL0AD*L/15 . ) 

DO 46 1=2,6 

IM1=I-1 

DO 46 J=1 ,IM1 

46 KEL( I ,J ,EN)=KEL( J.I.EN), 

GENERATE THE LOCAL ELEMENT' MASS MATRICES AND STORE THEM IN THE 
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C 3-DIMENSIONAL MASS ARRAY MEL (SEE FIG. 8-11 FOR THE EQUATION 
C USED). EACH PLANE OF THE 3-DIM. ARRAY MEL CONTAINS ONE LOCAL 
C ELEMENT MASS MATRIX. 

F=GAMMA( EN )*L/420 . 

P=70 . *F 
P2=2 .*P 
0=156. *F 
S=22.*L*F 
T = 54 . *F 
U=4 . *L*L*F 
V=13.*L*F 
U=3 . *L*L*F 
MEL( 1 ,1 ,EN )=P2 
MEL( 1 ,4 ,EN )=P 
MEL( 2 ,2 ,EN )=0 
MEL( 2 ,3 ,EN )=S 
MEL( 2 ,5 ,EN )=T 
MEL( 2 ,6 ,EN )=-V 
MEL( 3 ,3 ,EN )=U 
MEL( 3 .5 ,EN )=V 
MEL( 3 ,6 ,EN )=-U 
MEL( 4 ,4 ,EN )=P2 
MEL( 5 , 5 ,EN )=0 
MEL( 5 ,6 ,EN )=-S 
MEL( 6 ,6 ,EN )=U 
DO 47 1=2,6 
IM1=I-1 
DO 47 J=1 ,IM1 

47 MEL( I ,J,EN)=MEL( J.I.EN) 

C GENERATE THE TRANSFORMATION MATRIX R AND ITS TRANSPOSE RT . 
R( 1 ,1 )=COSA 
R( 1 ,2)=SINA 
R( 2 , 1 )=-SINA 
R( 2 ,2 )=COSA 
R( 3 ,3 )=1 . 

R( 4 .4 )=COSA 
Rt4 ,5)=SINA 
R(5,4)=-SINA 
R( 5 ,5 )=COSA 
R( 6 ,6 )=1 . 

DO 48 1 = 1 ,3 
DO 48 J=1 ,3 
RT( I , J )=R( J , I ) 

48 RT(I+3,J+3)=R( J+3,1+3) 

C DETERMINE THE ELEMENT STIFFNESS MATRICES IN THE GLOBAL 
C COORDINATE SYSTEM ( EQ . 8-93B ) AND STORE THEM IN THE 3-DIM. 

C STIFFNESS ARRAY KEG. EACH PLANE OF THE 3-DIM. ARRAY CONTAINS 
C ONE GLOBAL ELEMENT STIFFNESS MATRIX. 

DO 95 1 = 1 ,6 
DO 95 J=1 ,6 
TK( I , J )=0 .0 
DO 95 K=1 ,6 

95 TK( I , J )=TK( I ,J )+KEL( I ,K ,EN )*R( K , J ) 

DO 96 1=1,6 
DO 96 J=1 ,6 
KEG( I , J ,EN )=0 .0 
nn qa k = i .a 
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96 KEG( I , J ,EN )=KEG( I , J ,EN )+RT( I ,K )*TK( K , J ) 

C DETERMINE THE ELEMENT MASS MATRICES IN THE GLOBAL SYSTEM ( EO . 

C 8-93A) AND STORE THEM IN THE 3-DIM. MASS ARRAY MEG. EACH PLANE 
C OF THE 3-DIM. ARRAY CONTAINS ONE GLOBAL ELEMENT MASS MATRIX. 

DO 97 1*1 ,6 
DO 97 J = 1 ,6 
TM( I ,J)*0.0 
DO 97 K = 1 ,6 

97 TM( I . J )=TM( I . J ) + MEL( I # K , EN )*R( K f J ) 

DO 98 1 = 1 , 6 

DO 98 J*1 ,6 
MEG( I , J , EN )=0 . 0 
DO 98 K = 1 ,6 

98 MEG( I , J , EN )=MEG( I , J , EN )+RT( I ,K )*TM( K , J ) 

100 CONTINUE 

C GENERATE NULL MATRICES SK AND SM WHICH WILL BECOME THE 
C SYSTEM STIFFNESS AND MASS MATRICES , RESPECTIVELY. 

N=NJTS*3 
DO 49 1 = 1 ,N 
DO 49 J = 1 , N 
SK( I ,J )=0 . 

49 SM( I ,J )=0 . 

C ASSEMBLE THE STIFFNESS AND MASS MATRICES. 

DO 51 1 = 1 ,NUMEL 
DO 50 J = 1 ,2 
DO 50 M= 1 ,3 
J1=J*3-M+1 

50 SUB( J1 )=3*JNM( I ,J )-M+l 
DO 51 B=1 ,6 

DO 51 Z = 1 ,6 
ROWSUB=SUB( B ) 

COLSUB=SUB( Z ) 

SK( ROWSUB ,COLSUB )=SK( ROWSUB , COLSUB )+KEG( B ,Z , I ) 

51 SM( ROWSUB , COLSUB )=SM( ROWSUB t COLSUB )+MEG( B f Z , I ) 

C CALCULATE THE NUMBER OF DEGREES OF FREEDOM AND REMOVE ROWS AND 
C COLUMNS FROM THE SYSTEM STIFFNESS AND MASS MATRICES. 

NF=N-NB 

IF( NB .EQ. 0 )GQ TO 69 

NA=1 

KL=N-1 

62 JC=1 

63 IF( JC .EQ. CFIX(NA))GO TO 64 
JC=JC+1 

IF(JC .EQ. N )G0 TO 68 
GO TO 63 

64 DO 65 1=1 ,N 
DO 65 J*JC , KL 

SK( I # J )=SK( I , J + l ) 

SM( I ,J)=SM( I ,J + 1 ) 

65 CONTINUE 

DO 66 J*1 ,N 
DO 66 1 = JC ,KL 
SK( I , J )=SK( 1 + 1 f J) 

SM( I , J )=SM( 1 + 1 ,J) 

66 CONTINUE 

IF( NA .EQ. NB )G0 TO 68 
NA=NA+1 
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DO 67 I*NA,NB 

67 CF IX( I )*CFIX( I )-l 
GO TO 62 

68 CONTINUE 

C ASSIGN REDUCED STIFFNESS AND MASS MATRIX ELEMENTS TO ARRAY 

C NAMES RSK AND RSM , RESPECTIVELY. 

69 DO 70 1*1 ,NF 
DO 70 J*1 ,NF 
RSK( I , J )=SK( I ,J) 

70 RSM( I , J )*SM( I , J ) 

C WRITE OUT THE REDUCED STIFFNESS AND MASS MATRICES OBTAINED 

C FROM THE BOUNDARY CONDITIONS. 

WRITEC 2,71) 

71 FORMATC /,’ THE REDUCED SYSTEM STIFFNESS MATRIX IS:’,/) 
WRITEC 2,72) ( ( RSK( I ,J ),J*1 ,NF),I = 1 ,NF ) 

72 FORMAT( ’ ’ ,6EU .4/ ) 

WRITEC 2 ,73 ) 

73 FORMATC / , ’ THE REDUCED SYSTEM MASS MATRIX IS:’,/) 

WRITEC 2,74) (CRSM(I,J),J*1,NF),I«1,NF) 

74 FORMATC * ’.6E11.4/) 

WRITEC 2 ,200 ) 

200 FORMATC //,’’) 

C CALL SUBPROGRAM JCBI TO CALCULATE FREQUENCIES AND MODE SHAPES 
CALL JCBIC NF ,RSK ,RSM ) 

STOP 

END 

C LIBRARY. FOR 

C SUBROUTINES JCBI, DECOMP, MATINV , MATMPY , AND SEARCH AS REQU1 

C BY FINITEL .FOR AND TRUSS .FOR 
SUBROUTINE JCBIC N.K.M) 

IMPLICIT REAL*8 CA-H.O-Z) 

REAL*8 k ,m,l ,lt ,LINV,LINVTR,RT , a, omega, prod, av ,diff ,RAD 

REAL*8 COSINE, SINE, Q.PR0D1 

DIMENSION K( 27,27 ) ,RT( 27,27) ,A( 27,27 ) 

DIMENSION OMEGAC 27 ) ,MC 27 ,27 ) ,L( 27 ,27 ) ,LT( 27,27) 

DIMENSION LINVC 27 ,27 ) .LINVTRC 27.27 ) .PRODC 27,27) 

CALL DECOMPC M ,N ,L ,LT ) 

CALL MATINVCL.LINV.N) 

DO 204 1*1 ,N 
DO 204 J=1 ,N 

204 LINVTRC I , J )“LINVC J ,1 ) 

CALL MATMPYC N ,K .LINVTR ,PROD ) 

CALL MATMPYC N.LINV, PROD, A) 

DO 14 1*1, N 
DO 13 J*1,N 
RT( I , J )*0 .0 

13 CONTINUE 
RT( I ,1 )*1 .0 

14 CONTINUE 
NSWEEP-0 

15 NRSKIP-0 
NMIN1=N-1 

DO 25 1*1 ,NMIN1 

IP1*I+1 

DO 24 J*IP1,N 

AV=0 ,5*( A( I , J )+A( J ,1 ) ) 

nTFF.Ar T Tl-4n..l) 
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RAD=DSQRT( DIFF*DIFF+4 .*AV*AV ) 

ifcrad.eq.o.o )go to 20 

IF( DIFF .LT .0 .0 )G0 TO 18 

Go' T0 B 17* ( 2 ‘ 2 ^ -OABSC A( 1 ,1 ) )+100 .*DABS( Av ) )G0 TO 1G 

It co4?N^DSQRf( J c > R<;o?i?Fr?Sr2 J ifR2o 1 >° 0 ' 0 ' D ' ,B£< Av > >eo T0 20 

SINE=AV/( RAD*COSINE ) " 

GO TO 19 

18 ?i?fu D , s 2 RT( ( )/( 2 . 0*RAD ) ) 

IF(AV.LT.O.O)SINE=-SINE 
COSINE=AV/( RAD*SINE ) 

19 dIsSb* SINE > STaTEMENT ° R >«NA L program 

IF( DBS .GT .1 . OE-16 )G0 TO 21 

20 NRSKIP-NRSKIP+1 
GO TO 24 

21 DO 22 Ll = l , N 
Q=A( I ,L1 ) 

A( I ,L1 )=C0SINE*0+SINE*A( J ,L1 ) 

AC J ,L1 )=-SINE*0+C0SINE*A( J ,L1 ) 

22 CONTINUE 1 

DO 23 Ll = l ,N 

Q=A( LI .1 ) 

A( LI ,1 )=C0SINE*0+SINE*AC LI ,J ) 

A( LI , J )=-SINE*Q+COSlNE*A( LI . J ) 

Q=RT(L1,I) 

RT( LI , I )=COSINE*Q+SlNE*RTC LI , J ) 

RT( LI , J )=-SINE*Q+COSINE*RT( LI . J ) 

23 CONTINUE * 

24 CONTINUE 

25 CONTINUE 

KEEP A TALLY OF THE NUMBER OF SWEEPS 
NSWEEP=NSWEEP+1 
IF( NSUEEP .GT . 100 )G0 TO 33 
WRITEC 2,26 )NRSKIP .NSWEEP 

26 FORMATC • ’,5X, ’THERE WERE ’,12 

ROTATIONS SKIPPED ON SWEEP NUMBER * 12) 

IF( NRSKIP .LT ,N*( N-l )/2 )G0 TO 15 ' * 

PR0D1-0 .0 
DO 27 J=l ,N 

PR0D1=PR0D1+RT( J , 1 )*RT( J ,N ) 

27 CONTINUE 
WRITEC 2, 28) 

28 )PR00i ' THE 8CAL '' R PR0DU1:T 0F THE FIRST AND LAST ’> 
2 %F19 M 17/) ’ ,5X ' 'EIGENVECTORS OF THE TRANSFORMED MATRIX IS 

CALL MATMPYCN.LINVTR.RT .PROD) 

DO 30 1=1 ,N 
DO 30 J=1 ,N 
30 RT(I,J)«PR0DCl,J) 

DO 42 J«1 ,N 
SUM=0,0 


DO 31 1=1 ,N 

31 SUM-SUM+DABSCRTC I ,J)) 
AV-SUM/N 
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QUOT=DABSC RT( 1 ,J))/AV 

IF( QUOT .LT .0 .000001 )G0 TO 40 

DO 32 1=2, N 

32 RT( 1 , J )=RT ( I ,J )/RT( 1 , J ) 

RT( 1 , J )*1 .000 

GO TO 42 

40 CALL SEARCHC RT , J , 1 1 ,N ) 

BIG=RT( II ,J) 

DO 41 I«l , N 

41 RT( I ,J )=RT( I ,J )/BIG 

42 CONTINUE 

DO 110 1=1 ,N 

IF(A(I,I ).LE.0.0)G0 TO 43 
OMEGAC I )=DSQRT( A( I , I ) ) 

GO TO 110 

43 OMEGAC I )=0 .0 

110 CONTINUE 

33 WRITEC 2 ,34 )NSWEEP 

'34 FORMAT(/,’ ’,5X, 'THERE WERE ’,13,’ SWEEPS PERFORMED.'. 
*/,5X,' THE EIGENVALUES AND EIGENVECTORS FOLLOW:’) 

DO 39 JJ = 1 , N 
J =N- J J + 1 

WRITE( 2 ,35 )JJ ,A( J , J ) 

35 FORMAT( / , ’ ’,5X, 'LAMBDA (',12,') = ’.F20.4) 

WRITEC 2,111 )JJ ,OMEGA( J ) 

111 FORMAT( ’ ’ ,5X, 'OMEGA( ’ ,12, ' ) « ’.F20.4.* RAD/S’) 

WRITEC 2, 36) 

36 FORMAT( / , ’ ’ ,5X , 'THE ASSOCIATED EIGENVECTOR IS:*) 

DO 37 1=1, N 

37 WRITEC 2 ,38 )RT( I , J ) 

38 FORMATC ’ ’,5X,D17.10) 

39 CONTINUE 
RETURN 
END 


SUBROUTINE DECOMPC A.N.L ,LT) 
IMPLICIT REAL*8( A-H.O-Z) 
DIMENSION AC 27, 27) 

REAL*8 LC27,27),LTC27,27) 

DO 9 3*1, N 

IFCJ.EQ.DGO TO 7 

JM1=J-1 

DO 6 I=J,N 

IFC I . NE . J )G0 TO 4 

SUM=0 .0 

DO 3 K=1 , JM1 

3 SUM=SUM+LC I ,K)*L( J,K) 

LC J . J )=DSQRTC AC J . J )-SUM ) 

GO TO 6 

4 SUM=0.0 

DO 5 K*1 , JM1 

5 SUM=SUM+LCI,K)«LC J.K) 

L( I , J )*C AC I , J )-SUM )/LC'J , J ) 

6 CONTINUE 
GO TO 9 
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c 


c 


7 

8 
9 


L( 1 ,1 )=DSQRT( A( 1 . 

DO 8 1=2, N 

<-( I .1 )*A( I .1 )/|_( 1 

CONTINUE 


)) 

1 ) 


FIL D0 I ll 2 j-2 r UES ° F MATRIX L 

JM1 = J- 1 
DO 11 I = i ,JM1 

ASS 00 N 12 S ^UN TO 0PPER TRI ' wsu '-«< "OTBIX 
DO 12 J=1 In 
12 LT(I.J)«L(J,I) 

RETURN 

ENn 


LT 


C 

C 

C 


C 


C 


SUBROUTINE MATINV( B , A ,N ) 

MATRIX INVERSION USING GAUSS-JOROAN REDUCTION AND PARTIAL 
PIVOTING. MATRIX B IS THE MATRIX TO BE INVERTED AND A IS 
THE INVERTED MATRIX. 

IMPLICIT REAL*8( A-H.O-Z) 

DIMENSION B( 27 ,27 ) ,A( 27 .27 ) , INTER( 27 ,2 ) 

DO 2 1 = 1, N 


DO 2 J=1 , N 
2 A( I , J )=B( I , J ) 

CYCLE PIVOT ROW NUMBER FROM 1 TO N 
DO 12 K=1 ,N 
JJ=K 


IF(K.EO.N)GO TO 6 
KP1*K+1 

BIG=DABS( A( K ,K ) ) 

SEARCH FOR LARGEST PIVOT ELEMENT 
DO 5 I=KP1 ,N 
AB=DABS( A( I ,K ) ) 

IF( BIG-AB )4 , 5 , 5 
4 BIG=AB 
JJ*I 
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5 CONTINUE 

C MAKE DECISION ON NECESSITY OF ROW INTERCHANGE AND 

C STORE THE NUMBER OF THE TWO ROUS INTERCHANGED DURING KTH 

C REDUCTION. IF NO INTERCHANGE, BOTH NUMBERS STORED EQUAL K 

6 INTER( K , 1 )«K 
INTER( K ,2 )=JJ 
IF( JJ-K )7 ,9 , 7 

7 DO 8 J = 1 ,N 
TEMP=A( JJ ,J ) 

A( JJ , J )=A( K , J ) 

8 A(K,J)=TEMP 

C CALCULATE ELEMENTS OF REDUCED MATRIX 

C FIRST CALCULATE NEU ELEMENTS OF PIVOT ROU 

9 DO 10 J = 1 ,N 
IF(J.EO.K)GO TO 10 
A(K.J)=A(K,J)/A(K,K) 

10 CONTINUE 

C CALCULATE ELEMENT REPLACING PIVOT ELEMENT 
A( K , K )=1 . /A( K , K ) 

C CALCULATE NEW ELEMENTS NOT IN PIVOT ROU OR COLUMN 
DO 11 1 = 1 ,N 
IF( I .EQ.K)GO TO 11 
DO 110 J=1 ,N 
IF(J.EQ.K)GO TO 110 
A( I , J )=A( I , J )-A( K , J )*A( I ,K ) 

110 CONTINUE 

11 CONTINUE 

C CALCULATE NEU ELEMENTS FOR PIVOT COLUMN— EXCEPT PIVOT ELEMENT 
DO 120 1 = 1 , N 
IF(I.EQ.K)GO TO 120 
A( I , K )=— A( I ,K )*A( K ,K ) 

120 CONTINUE 

12 CONTINUE 

C REARRANGE COLUMNS OF FINAL MATRIX OBTAINED 
DO 13 L=1 , N 
K=N-L+1 

KROW=INTER( K , 1 ) 

IROU=INTER( K ,2 ) 

IF( KROU .EQ . IROU )G0 TO 13 
DO 130 1=1 ,N 
TEMP=A( I .IRON) 

A( I , IROU )=AC I .KROU 
A( I ,KROU)=TEMP 
130 CONTINUE 

13 CONTINUE 
RETURN 
END 

C 

SUBROUTINE MATMPY( N , A ,B ,C ) 

IMPLICIT REAL*8( A-H.O-Z) 

C C IS THE PRODUCT MATRIX OF A AND B 

DIMENSION A( 27 ,27 ) ,B( 27 ,27 ) ,C( 27 ,27 ) 

DO 2 1*1 ,N 
DO 2 J=1,N 
C(I ,J)*0.0 
DO 2 K=1 ,N 

2 C( I , J )=C( I , J )+A( I ,K )*B( K ,J ) 

RETURN 

END 
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SUBROUTINE SEARCH( RT , J , 1 1 ,N ) 

C THIS SUBROUTINE SEARCHES THE JTH COLUMN 
C FOR THE LARGEST EIGENVECTOR COMPONENT. 

C ASSIGNED TO THE NAME II. 

IMPLICIT REAL*8( A-H.O-Z) 

DIMENSION RT( 27 , 27 ) 

11 = 1 

BIG=DABS( RT( 1 ,J)) 

DO 3 1=2, N 
AB=DABS( RT( I , J ) ) 

IF( BIG-AB )2 ,3 ,3 

2 BIG=AB 
II = I 

3 CONTINUE 
RETURN 
END 


OF THE MATRIX RT 
ITS ROU NUMBER IS 
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NLFINITE .FOR; FULL KE+KG MATRICES .REVISED 10-24-90 
THE AREA ARRAY A IS: 

A( 1 ) « 0 .4800000E+02 

THE ELASTICITY ARRAY E IS; 

E( 1 ) * 0 .3000000E+08 

THE MOMENT OF INERTIA ARRAY IA IS; 

IA( 1 ) = 0 . 1000000E+04 

THE GAMMA ARRAY IS; 

GAMMA(l) = 0.3525000E-01 

THE AXIAL PRETENSION LOAD IS; 

0 . 

THE JOINT-NUMBER MATRIX IS; 

1 2 

THE JOINT COORDINATES ARE; 

X( 1 )= O.OOOOOOOE+OO Y( 1 )* O.OOOOOOOE+OO 
X( 2 )= 0 . 1000000E+03 Y( 2 )= O.OOOOOOOE+OO 

THE REDUCED SYSTEM STIFFNESS MATRIX IS: 

0 . 1440E+08 O.OOOOE+OO 0 .OOOOE+OO-O . 1440E+08 O.OOOOE+OO O.OOOOE+OO 

O.OOOOE+OO 0 .3600E+06 0.1800E+08 0 . OOOOE+OO-O .3600E+06 0.1800E+08 

O.OOOOE+OO 0.1800E+08 0.1200E+10 0 .OOOOE+OO-O . 1800E+08 0.6000E+09 

-0.1440E+08 O.OOOOE+OO O.OOOOE+OO 0.1440E+08 O.OOOOE+OO O.OOOOE+OO 

0 . OOOOE+OO-O . 3600E+06-0 . 1800E+08 O.OOOOE+OO 0 . 3600E+06-0 . 1800E+08 


O.OOOOE+OO 0 . 1800E+08 0.6000E+09 0 .OOOOE+OO-O . 1800E+08 0.1200E+10 
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THE REDUCED SYSTEM MASS MATRIX IS: 

0 . 1 1 75E+01 O.OOOOE+OO O.OOOOE+OO 0.5875E+00 O.OOOOE+OO O.OOOOE+OO 
0 . OOOOE+OO 0 . 1309E+01 0.1846E+02 O.OOOOE+OO 0 .4532E+00-0 . 1091E+02 
0 . OOOOE+OO 0 . 1846E+02 0.3357E+03 O.OOOOE+OO 0 . 1091E+02-0 . 251 8E+03 
0 . S875E+00 O.OOOOE+OO O.OOOOE+OO 0.1175E+01 O.OOOOE+OO O.OOOOE+OO 
O.OOOOE+OO 0 . 4532E+00 0.1091E + 02 O.OOOOE+OO 0 . 1309E+01-0 . 1846E+02 
0 .OOOOE+OO-O .1091 E+02-0 .251 8E+03 0 .OOOOE+OO-O .1846E+02 0.3357E+03 


THERE UIERE 5 ROTATIONS SKIPPED ON SWEEP NUMBER 1 
THERE WERE 9 ROTATIONS SKIPPED ON SWEEP NUMBER 2 
THERE WERE 9 ROTATIONS SKIPPED ON SWEEP NUMBER 3 
THERE WERE 14 ROTATIONS SKIPPED ON SWEEP NUMBER 4 


THERE WERE 15 ROTATIONS SKIPPED ON SWEEP NUMBER 5 
THE SCALAR PRODUCT OF THE FIRST AND LAST 

EIGENVECTORS OF THE TRANSFORMED MATRIX IS 0.00000000000000000 


THERE WERE 5 SWEEPS PERFORMED. 

THE EIGENVALUES AND EIGENVECTORS FOLLOW: 

LAMBDA ( 1) ■ 0.0000 

OMEGA( 1) * 0.0000 RAD/S 

THE ASSOCIATED EIGENVECTOR IS: 

0 .1000000000D+01 
O.OOOOOOOOOOD+OO 
O.OOOOOOOOOOD+OO 
0.1000000000D+01 
O.OOOOOOOOOOD+OO 
O.OOOOOOOOOOD+OO 

LAMBDA ( 2) * 0.0000 

OMEGA( 2) « 0.0000 RAD/S 

THE ASSOCIATED EIGENVECTOR IS: 
O.OOOOOOOOOOD+OO 
0 .9334669755D+00 
0.6653302446D-03 
0 .OOOOOOOOOOD+OO 
0.1000000000D+01 
0.6653302446D-03 



LAMBDA ( 3) 

OMEGA( 3) = 


0.0000 

0.0000 RAD/S 


THE ASSOCIATED EIGENVECTOR IS: 

0 .OOOOOOOOOOD+OO 
0.1000000000D+01 
-0.1977319320D-01 
0. OOOOOOOOOOD+OO 
-0 .97731 93204 D+OO 
-0 . 19773 19320D- 01 

LAMBDA ( 4) = 6127659.5745 

OMEGA( 4) * 2475.4110 RAD/S 

THE ASSOCIATED EIGENVECTOR IS: 

0 .OOOOOOOOOOD + OO 
0.1000000000D+01 
-0.6000000000D-01 
0 .OOOOOOOOOOD+OO 
0.1000000000D+01 
0 .60000000000-01 


LAMBDA ( 5) = 49021276.5957 

0MEGA( 5) = 7001.5196 RAD/S 

THE ASSOCIATED EIGENVECTOR IS: 
0.1000000000D+01 
0. OOOOOOOOOOD+OO 
0. OOOOOOOOOOD+OO 
-0.1000000000D+01 
0 .OOOOOOOOOOD+OO 
0. OOOOOOOOOOD + OO 


LAMBDA ( 6) = 7 1489361.7021 

OMEGA( 6) = 8455.1382 RAD/S 


THE ASSOCIATED EIGENVECTOR IS: 
0. OOOOOOOOOOD+OO 
0.1000000000D+01 
-0 . 1200000000D+00 
0 .OOOOOOOOOOD+OO 
-0.1000000000D+01 
-0 .1200000000D+00 
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APPENDIX D 

3 -Node Beam Derivation of [Kg] 
Consider the 3 -node beam shown in Figure 23 . 

2 

u(x) = aQ + a^x + a2X 
u' = a]_ + 2a2X 
u" = 2a 2 

(u') 2 = a^ 2 + 4aia2X + 4a2 2 x 2 

2 3 4 

v(x) = bg + b^x + b2X + b3X + b4X 
v ' = bj + 2b2X + 3 b 3 X 2 + 4 b 4 X 3 
v" = 2b2 +6b3X +12b4X 2 

(v') 2 = 16x^b4 2 + 24 x^b 3 b 4 + 16x 4 b2b4 + 9x 4 b3 2 + 8x 3 b]_b4 
+ 12x^b2b3 + 6x 2 b^b3 + 4 x 2 b 2 2 + 4 xb^b 2 + b]_ 2 

uj = u(-L/ 2 ) = ag - a^L /2 + a2L 2 /4 
U2 = u( 0 ) = ao 

U3 = uL/ 2 ) = aQ + a^L /2 + a2L 2 /4 

v x = v(-L/2 ) = b 0 - biL /2 + b 2 L 2 /4 - b 3 L 3 /8 + b 4 L 4 /16 
v 2 = v(0) = b 0 

V3 = v(L/2 ) = bo + b^L /2 + b2L 2 /4 + b3L 3 /8 + b 4 L 4 /16 
0 ! = v'(-L/2 ) = bi - b 2 L + 3 b 3 L 2 /4 -b 4 L 3 /2 
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Figure 23 3-Node Beam Element 
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0 2 = v'(0) = b 1 

0 3 = v ' ( L/ 2 ) = bi + b 2 L + 3b 3 L 2 / 4 + b 4 L 3 /2 


Solve for and b^ 


a 0 = u 2 
bo = v 2 

b ! = ®2 


2 

u l - u 2 = -a^L/2 + a 2 L /4 

2 

u 3 ~ u 2 = a l L /2 + a 2 L /4 

2 

U 1 " 2u 2 + u 3 = a 2 L /2 


a 2 


2u^ - 4 u 2 + 2u 3 



U 1 - u 3 = _a l L 


u 3 - u x 

a l = 

L 


3 

v 3 “ V 1 " = ^3 l /4 


4V3 4v^ 402 



3 

v 3 " V 1 = b^L + b3L / 4 

V 1 + v 3 " 2 v 2 = 2bo + b2L^/2 + b4L^/8 - 2bg 




196 


®3 " ®1 = 2b2L + b4L^ 

4 vi/L + 4V3/L - 8V2/L = 2b2L + b4L^/2 
4 vi/L + 4V3/L - 8V2/L - 83 + 0i = -b4L^/2 



83 81 4 vi 4 v 3 8v 2 83 01 





197 



Integration yields 


L 

J[X] dx 
0 


L 2 4L 3 / 3 

0 0 

0 0 

0 0 


0 

0 

0 

0 

0 

0 

0 

0 

L 

L 2 

L 3 

L 4 

L 2 

4L 3 / 3 

3L 4 / 2 

8L 5 / 5 

L 3 

3L 4 / 2 

9L 5 /5 

2L 6 

L 4 

8L 5 /5 

2L 6 

16L 7 /7 


0 


0 
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a l a 2 b l b 2 b 3 b 4 


[L!] T 

U 1 V 1 ®1 u 2 v 2 ®2 u 3 v 3 83 J "1/L 2/L 2 00 0 0 

000 4/L 2 -4/L 3 - 8 /L 4 

000 1/2L 0 -2/L 3 

0 -4/L 2 000 0 

000 - 8 /L 2 0 16/L 4 

0 010 -4/L 2 0 

1/L 2/L 2 000 0 

000 4/L 2 4/L 3 - 8 /L 4 

000 -1/2L 0 2/L 3 



7 

3L 

0 


0 


-20 

3L 


0 


0 

13 

3L 

0 


0 


0 

0 

-20 

3L 

0 

0 

13 

3L 

0 

0 

18272 

3469 

0 

-22096 

304 

0 

3824 

-3469 

105L 

105 

105L 

5 

105L 

105 

3469 

659L 

0 

-4208 

23L 

0 

739 

-659L 

105 

105 

105 

2 

105 

105 

0 

0 

64 

3L 

0 

0 

-44 

3L 

0 

0 

-22096 

-4208 

0 

27392 

-72 

0 

-5296 

4208 

105L 

105 

105L 

105L 

105 

304 

23L 

0 

-72 

109L 

0 

56 

-23L 

5 

2 

5 

5 

2 

0 

0 

-44 

3L 

0 

0 

31 

3L 

0 

0 

3824 

739 

0 

-5296 

56 

0 

147 

-739 

105L 

105 

105L 

5 

105 

105 

-3469 

-659L 

0 

4208 

-23L 

0 

-739 

659L 

105 

105 

105 

2 

105 

105 



[Kg] 

3 -NODE 

ELEMENT 
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3 NODE RIGID BODY ROTATION VECTOR 


Exact 


Expanded to 


Factor out 


2 Terms 


BL 


L/2 (1-cos 2fl) 


L(fl 2 -B 4 /3) 


B-B 3 / 3 

-L/2 sin 2B 


-L ( B-2B 3 / 3 ) 


-1 + 2B 2 / 3 

2fl 


2B 


2/L 

0 


0 


0 

0 


0 


0 

2B 


2B 


2/L 

-L/2 (1-cos 26) 


-L(B 2 -fl 4 /3) 


-fl + fl 3 / 3 

L/2 sin 2B 


L(B-2B 3 /3) 


1 - 2B 2 /3 

2B 


2B 


2B 


91.733 B 3 - 16 B 
L( 17 . 333 B 3 - 3 B 
8 B 2 

16 B - 106.666 B 3 
L(33.0666 B 3 - 6 B 
- 6 B 2 
14.9333 fl 3 
L ( 3 fl - 17.3333 fl 3 ) 


RFORCES using expanded 
rigid body rotation vector 
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cos 2B - 1 

121.6 B - 68.6 sin 2B 
23 LB - 13L sin 2B 
4-4 cos 2B 

80 sin 2B - 144 B 

43.6 LB - 24.8 L sin 2B 
3 cos 2B - 3 

22.4 B - 11.2 sin 2B 
13 L sin 2B - 23 L B 


RFORCES using exact 
rigid body rotation 
vector 
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Abstract 

In order to be cost-effective, space structures must 
be extremely light-weight, and subsequently, very flexible 
structures. The power system for Space Station Freedom is 
such a structure. Each array consists of a deployable truss 
mast and a split "blanket" of photo-voltaic solar 
collectors. The solar arrays are deployed in orbit, and the 
blanket is stretched into position as the mast is extended. 
Geometric stiffness due to the preload make this an 
interesting non-linear problem. 

The space station will be subjected to various 
dynamic loads, during shuttle docking, solar tracking, 
attitude adjustment, etc. Accurate prediction of the 
natural frequencies and mode shapes of the space station 
components, including the solar arrays, is critical for 
determining the structural adequacy of the components, and 
for designing a dynamic controls system. 

This paper chronicles the process used in developing 
and verifying the finite element dynamic model of the photo- 
voltaic arrays. Various problems were identified in the 
investigation, such as grounding effects due to geometric 
stiffness, large displacement effects, and pseudo-stiffness 
(grounding) due to lack of required rigid body modes. 
Various analysis techniques, such as development of rigorous 
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solutions using continuum mechanics, finite element solution 
sequence altering, equivalent systems using a curvature 
basis, Craig-Bampton superelement approach, and modal 
ordering schemes were utilized. This paper emphasizes the 
grounding problems associated with the geometric stiffness. 


a 


Di 

d/dx, or 
d/dt, or 
E 


e a 

(F) 

F( x , t ) 

g 

i 

[K] 

[Ke] 

[Kg] 

L 


M 

dM 


m 

P 


P' 


Nomenclature 

factor defined by Eq.(13) 

arbitrary constants in Eq. (10) 

differential operator with respect to position 

differential operator with respect to time 

modulus of elasticity 

axial strain 

input force vector at the beginning of a step 

applied transverse force 

factor defined by Eq.(14) 

moment of inertia 

stiffness matrix 

elastic stiffness matrix 

geometric stiffness matrix 

1 ength 

moment 

change in moment 
mass per unit length 
axial force 

pseudo-force necessary for equilibrium 
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{R} 

T 

{u} 

u 

UA 

U B 

v 

V 
dV 

V 

dVol 

x 

y 

B 

6 

e 

e 


force vector, output force vector at the end of 
a step 

kinetic energy 

displacements at the node points 
longitudinal displacement 
strain energy due to axial load 
strain energy due to bending 
transverse displacement 
shear 

change in shear 

potential of the external loads 
change in volume 
axis defined by Figure 25 
axis defined by Figure 25 
1/2 the angle of rotation 
factor defined by Eq.(ll) 
factor defined by Eq.(12) 
angle of rotation 
stress 


Introduction 

NASA's Space Station Freedom consists of various 
modules supported by a space truss. Power for the space 
station will be provided by a deployable system of split 
blanket photo-voltaic arrays, which will have two degree of 
freedom rotational capabilities in order to track the sun 
during its orbit. The arrays are designed to be operated in 
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during its orbit. The arrays are designed to be operated in 
a zero-gravity environment. 

NASA Lewis Research Center, along with its 
contractors, have the responsibility for developing a 
verified finite element dynamics model of the solar arrays, 
which could be combined with the other space station 
substructures for both structural and dynamic control 
studies. The development of the model necessitated the use 
of unique procedures, and rigorous analytical checks. 

The procedure included the following: 

1. Development of an idealized model of the solar arrays, 
and derivation of a unique solution for the response 
frequencies for the idealized array cantilevered from 
the space truss, using equations developed from 
continuum mechanics .[ 1 ] 

2. Comparison of the frequencies from the MSC/NASTRAN 
finite element dynamic model of the idealized array 
with the rigorous solution from continuum mechanics .[ 2] 

3. Refinement of the finite element mesh. 

4. Rigid body mode checks of the finite element models. 

5. Various parameter studies involving the amount of 
tension in the blanket, rigidity of the blanket tip 
beam, type of elements used, etc.. 
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6. Craig-Bampton approach for appending rigid body modes 
to substructures (superelements) [3], 

7. Modal ordering schemes for identifying "important" 
modes . 

8. Study of grounding effects due to lack of rigid body 
mode capabi 1 ities . [ 4 ] 

A detailed summary of the project was presented [5], 
It should be noted that this study is ongoing at the present 
time. This paper will be restricted to the grounding 
problems associated with the geometric stiffness due to 
blanket pre-load. 


Grounding 

The space station solar arrays were modeled 
utilizing MSC/NASTRAN. As a routine check, the stiffness 
matrices generated by the model were multiplied by a matrix 
of rigid body modes, and large pseudo-forces were developed 
(grounding). The cause of this "grounding" phenomenum was 
examined. 

Finite element solves non-linear problems of the 

form 

[[Ke] + [Kg]] * {u} = (R> - {F} 

where [Ke] is the elastic stiffness matrix, and [Kg] is the 
geometric, or initial stress stiffness matrix. 
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[Kg] is a function of the pre-load. Thus, it equals 
zero for a linear problem. [Ke] possesses the required 
rigid body modes. However, [Kg] lacks the capacity for 
rigid body rotation. Hence, an erroneous stiffening, or 
"grounding", occurs when a pre-loaded beam deforms. 

The traditional, or consistent geometric stiffness 
matrix, developed by Martin [6] and others, is 


6/5L 

1/10 

-6/5L 

1/10 

1/10 

2L/15 

-1/10 

-L/30 

-6/5L 

-1/10 

6/5L 

-1/10 

1/10 

-L/30 

-1/10 

2L/15 


This matrix does not possess rigid body rotation 
capabilities. Various refinements to the geometric 
stiffness have been developed which contain higher order 
terms [6,7,8] . However, none of these possess all the 

rigid body modes. Bosela [4] developed a modified [Kg] with 
complete rigid body modes when used with an exact rigid body 
rotation matrix, but [Kg] lost some of its rigid body 
capabi 1 i ties . 

Closer examination of the traditional formulation of 
[Kg] indicated that there is a load imbalance in the 
representation, and that pseudo-forces occur to maintain 
equilibrium (Figure 24). 
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the lack of rigid body rotation capabilities for [Kg] is not 
a problem, because the energy representation is correct. It 
can be shown that it is correct to B* terms, but error does 
occur, as a function of fl . For large rigid body rotation, 
as will occur with the solar arrays, this is significant. 

It should be noted that as long as the pre-load P is 
assumed to remain horizontal during rotation, work will be 
done by the force. Thus, true rigid body rotation cannot 
occur. In order for the strain energy to equal zero, the 
force P must change its orientation as the beam rotates ( 
ie. a follower force). 

Rigorous Solution Of Pre-Loaded Beam 

Suppose we have an axially loaded beam in space 
subjected to a time varying transverse loading (Figure 25). 

L 

The kinetic energy is f m (v')* 

T = dx (1) 

0 

The strain energy due to bending is 

r e i 

u B = — (v")‘ dx (2) 
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y,v 



dx 



Figure 25 


Beam in Tension and Differential Element 
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The strain energy due to axial load is 

U A = - | o e a dVol (3) 

Letting dVol = dA dx and applying non-linear elasticity 
yields 



(du/dx)’ 


+ du/dx(dv/dx) * 


+ l/4(dv/dx) |dx 


(4) 


Neglecting axial displacement and higher order terms yields 



(5) 


The potential of the external loads is 


V = - 


F( x , t ) v dx + V D v(0 , t) + M 0 v ' ( 0 , t ) 
- V L v(L , t ) - M l v ' (L , t ) 


( 6 ) 


Applying Hamilton's principle, and performing the 
variation, yields 
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t 2 L 



tl o 


£eIv"6(v")+Pv' 6( v' )-mv6(v)-F(x, t)6(v)J dx 

+V o 6v(0,t) +M 0 6 v ' (0,t)-V L 6v(L,t) -M L 6 v ' ( L , t ) j dt = 0. 

(7) 


Integrating by parts yields the differential equation 


d*/dx* (EId*v/dx* ) - P d'v/dx’ + m d*v/dt* = F(x,t) , ( 8 ) 


which agrees with Clough in reference [10], after a sign 
change required to express the axial force in tension 
instead of compression. This is also in agreement with 
Shaker in Reference [11]. 

For a beam in space, the moment and shear at the end 
points must equal zero. Thus, the boundary conditions are 

EIv"(0,t)=EIv"(L,t)=v'"(0,t)-P v’ (0,t)=v , "(L,t)-Pv' (L,t)=0 

El El (9) 


Choose a solution of the form 
v(x) = Disin(Sx) + D 2 Cos( 6 x) + D 3 sinh(ex) + D 4 COsh(ex). 

( 10 ) 


where 


& 


(a 4 +g 4 /4) 1 / 2 -g'/2] 


( 11 ) 
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e = [(a 4 +g 4 /4) 1 / 2 +g’/2] 

a 4 = mw* /El 
g* = P/EI 


( 12 ) 

(13) 

(14) 


Applying the boundary conditions at x = 0, and after much 
mathematical manipulation, yields 


v(x) 


D3 


6 sin 8 x + 
e 


sinhex 


+ 


D 4 


e * cosSx 
8 ' 


+ 


coshex 


(15) 


Applying the boundary conditions at x=L, and after more 
mathematical manipulations, yields 


D 3 ^ 6 ^cosheL - 8 3 cos 8 Lj + D 4 £ 3 sin 8 L + 8 3 sinh£Lj . (16) 


Expressing Eq.(15) and Eq.(16) into matrix form, setting the 
determinant equal to zero, and after more mathematical 
manipulations, the following characteristic equation is 
obtained 

±2a 6 (cosheLcos5L-l ) + ( £^- 6 ^)sinh£Lsin 8 L = 0 . (17) 

Using Eq.(13), this can be expressed as 

±w 3 (m/EI) 3 / 2 (cosh£Lcos 6 L-l) + ( £ 6 - 8 6 )sinh£Lsin 8 L = 0 . (18) 

By observation, when w=0, a=0, and 8=0. Letting 


sin( 0)=0 yields 
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w^ (m/EI ) ^ ^ ( coshsLcosSL-l ) = 0 


( 19 ) 


3 

The w term indicates that there must be three zero roots of 
"w" , which suggests the three required rigid body modes. 


Conclusion 

Lack of complete rigid body mode capabilities is 
inherent in the physical representation of the pre- tensioned 
beam problem currently used to formulate the geometric 
stiffness matrix. This lack of complete rigid body mode 
capabilities invalidates the rigid body mode check for non- 
linear problems, and adversely impacts the use of 
traditional finite element techniques to predict dynamic 
response of pre-loaded structures unless the missing rigid 
body modes are somehow apppended on to the structure, such 
as by the Craig-Bampton technique. 

The rigorous solution of the axial 1 y- 1 oaded beam 
with free/free boundary conditions developed in this paper 
may lend itself to the development of a new geometric 
stiffness matrix for a beam element with full rigid body 


capabi 1 ities . 
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APPENDIX F 

Diagonal izat ion/Partitioning Methodology 

Example 1 ( Structural Anal ysis . Third Edition, Ghali and 
Neville, Chapman and Hall Publishing Company, page 750.) 

Consider the beam in Figure 26. 


El 

[K] = 


1.6154 

-3.6923 

2.7692 


-3.6923 

10.1538 

-10.6154 


2.7692 

-10.6154 

18.4615 


[M] 


W 

g 


4 o 

o 1 

o o 


o 

0 

1 


{P} = P 0 [ 2, 1, 1 ] T sin fit 

[M] {’x } + [K] {X} = {0} 

| [K] - Q 2 [M] | = {0} 
yields 
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4W 0 2 Pq sin Qt 

L 


W o Pq sin Qt 

L 

W o ^ Pq sin Qt 

L 


77 


/ 


Figure 26 Example 1 
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, gEi 

Ql = 0.02588 — 3 D(x) = 

wL 13 


, gEi 

0 2 = 3.09908 — 3 D( 2 ) = 

wL -3 


, gEi 

1*3 = 25.89415 —3 D ( 3 ) = 

wL 


In matrix form, the Eigenvectors are 


[ 0 ] 


1.0 1.0 

0.5224 -6.3414 

0.1506 -4.5622 


We can use this transformation matrix 

A A 

[K] and [M] matrices. 


[K] = {0} T [K] {0} 

A 

[M] = {0} T [M] {0} 


1.0 

0.5225 

0.1506 

1.0 

-6.3414 

-4.5622 

1.0 

-13.1981 

19.2222 


1.0 

-13.1981 

19.2222 

to create diagonal 
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But, since 
Then, [K] 

CM] = 

[K] = Q 2 

In normal 

CM] {'n> + 
Where 


{’A} + [8] 


[K] {0} = Q 2 [M] {0}, 

A 

= a 2 [m] 

4.296 0 0 

0 65.027 0 

0 0 547.68 

4.296 0 

0 65.027 

0 0 


0 

0 

547.68 


coordinates, the equation of motion becomes 


[2] [M] {n> = {f} 


[0] 


Qi 2 0 0 

0 0 2 2 0 

0 0 Q 3 2 


and 


in = {0} T (p) 

= P 0 [ 2.673, -8.9036, 8.02412] T 


{n} = [M]" 1 {f} 
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or 


n l 


0.02588 0 0 


n l 


0.62206 

n 2 

gEI 

0 3.09908 0 


n 2 

P09 

-0.136921 


WL 3 



H 


n 3 


0 0 5.0886 


*3 


0.0146511 


Note that the equations are now un-coupl ed . 


Example 2 

Consider the beam in Figure 27. 
Let 

A = 48 in 2 
E = 30 x 10® psi 
I = 1000 in 4 
L = 100 in 

m = 0.03525 lb-sec 2 /in 2 


[Ke] = 


0 . 144xl0 8 
0 0 . 36xl0 6 

0 0.18xl0 8 

-0.144x10® 0 

0 -0.36xl0 6 

0 0.18x10® 


SYMMETRIC 

0.12X10 10 

0 0.144x10® 

-0.18x10® 0 0.36x10® 

0 . 6xl0 9 0 -0.18x10® 


.12x10 


10 
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PO 


sin fit Pq sin Qt 



Figure 27 Example 2 
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1.175 

0 1.309 SYMMETRIC 

0 18.46 335.7 

[M] = 

0.5875 0 0 1.175 

0 0.4532 10.91 0 1.309 

0 -10.91 -251.8 0 -18.46 335.7 


{P} = Pq [ 0 r 1, 0, 0, 1, 0 ] sin Qt 

2 

Ql = 0 

q 2 2 = 0 

Q 3 2 = 0 

Q 4 2 = 6,127,660 

Q 5 2 = 49,021,277 

Q 6 2 = 71,489,362 

1 0 0 0 1 0 

0 1 -50 1 0 1 

00 1 - 0.6 0 - 0.12 

C0] = 

1 0 0 0 -1 0 

0 1 50 1 0-1 


0 0 


1 


0.06 0 -0.12 
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[K] = [0] T [K] [0] = 



Example 3 

Now consider the beam with an axial load as shown in 
Figure 28. 

Let 

A = 48 in 2 
E = 30 x 10^ psi 
I = 1000 in 4 
L = 100 in 

m = 0.03525 lb-sec 2 /in 2 


T = 10,000,000 lbs 
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Pq sin Qt 


10,000,000 lbs 


L 


Pq sin Qt 


10,000,000 lbs 


Figure 28 Example 3 
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0 . 144xl0 8 

0 0.48x10® SYMMETRIC 

0 0.19x10® 0 . 1333xl0 10 

[ Ktan ] = 

-0.144x10® 0 0 0.144x10® 

0 -0 . 48x10® -0.19x10® 0 0.48xl0 6 

0 0.19x10® 0 . 5667xl0 9 0 -0.19x10® .1333xl0 10 

( p ) = Po [ o, 1, 0, 0, 1, 0 ] sin Qt 

1 0 0 0 1 

01-1 1 0 

0 0 0.018943 -0.06 0 

[03 = 

1 0 0 0 -1 

oil 10 

0 0 0.018943 0.06 0 

A 

[K] = {0} T [K] {<*} 

A _ 

[K] = 0 0 

0 0 
0 0 
0 0 
0 0 


0 0 0 0 

0 0 0 0 

4.0403 x 10 5 0 0 1226 

0 5.5174 x 10 6 0 0 

0 0 5.76 x 10 7 0 

3.8053 x 10 7 


0 

1 

-0.119554 

0 

-1 

-0.119554 


0 0 


1228 


0 


0 
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[ K ] tan contains a large erroneous term (4.0403 x 10 5 ) in 
the 3,3 position. 

Thus , 

: * 

[ R 3Tan * s not the correct diagonal matrix. 


The lack of rigid body rotation capability of the [Kg] 
matrix ultimately results in a large erroneous term in the 

A 

[K] matrix. 




